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I. 


INTRODUCTION 


A.  PROBLEM  CONTEXT 

Estimating  the  direction  of  arrival  (DOA)  of  a  signal  is 
a  problem  of  great  importance  in  the  operation  theater  of  a 
task  force  at  sea.  The  direction  of  arrival  (DOA)  vector  may 
be  used  to  distinguish  a  friendly  approaching  aircraft  from  a 
foe.  If  the  aircraft  is  outside  an  expected  bearing  corridor, 
then  it  is  generally  thought  to  be  an  enemy  aircraft.  This 
expectation  must  be  verified  by  other  check  points  during  the 
target  identification  and  classification  process.  Accurate 
determination  of  the  DOA  of  a  signal  is  a  very  important  issue 
in  the  operation  of  weapon  systems  and  should  be  completed 
with  extreme  care. 

The  main  goal  of  this  work  is  to  investigate  the  accuracy 
and  speed  of  high  resolution  DOA  techniques.  The  technique  of 
choice  will  allow  the  DOA  to  be  solved  accurately  in  real  time 
by  an  on  board  computer. 

Passive  linear  equispaced  phased  array  sensors  are  used  to 
generate  signals  for  the  DOA  determination.  Initially,  any 
signal  received  by  the  sensors  is  processed  to  estimate  signal 
direction  with  consecutive  samples  tracking  potentially  moving 
targets. 
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B .  THEORETICAL  BACKGROUND 

Consider  the  case  of  a  linear  equispaced  array  of  n 
sensors  receiving  signals  emitted  from  m  sources. 
Furthermore,  let  us  assume  that  the  signal  may  represented 
by  a  narrowband  signal  embedded  in  additive  Gaussian  white 
noise.  The  resulting  signal  received  at  the  array  is 

xc=s c+nc  ,  (1) 


where  s,  represents  all  narrowband  components  emitted  by  the 
source  (s)  and  is  the  additive  noise. 

The  output  at  the  q*  array  element  may  be  expressed  [Ref. 
1]  as 


A  r.r  2n  (g-l)  -d-sin  0, 

Aiexp[j-[wct - r - i+GJ]  +nCil 

i-1  A 


(2) 


for  1  s  g  s  n 


where  represents  the  i*  signal  arrival  angle,  wc  represents 
the  center  frequency  of  the  narrowband  sources,  d  represents 
the  array  element  spacing,  Aj  is  the  amplitude  of  each  incoming 
signal,  X  is  the  wavelength  related  to  the  traveling  wave 
movement,  represents  the  random  phase  of  each  incoming 
signal,  assumed  uniformly  distributed  over  [0,2 tt]  .  Finally, 
i\q  is  a  zero  mean  random  variable  representing  the  noise  at 
the  q*  array  element.  The  output  vector  at  time  t  is  a  n- 
dimensional  column  vector 


2 


(3) 


*£= 


'Vt.J 
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The  mode  vector  mg  is 

rriQ=  [1,  exp  ( -2x  j  cfsin  (0)  /X)  .  .  exp  ( -2 nj  d  (n-1)  sin(0)  /X)  ]  T. 

The  noise  vector  is 

nj.=  [nCil,  .  .  .  ,ntin]T.  (5) 


Equations  3-5  may  be  used  to  write  the  received  signal  vector 
as 


*£=£  A^exptj  [vet+*i]  ] 


1=1 


exp  ( -2it  j  dsin^)  /X) 


exp(-2x  j d  (n-1)  sin(0i)  /X) 


+nJ., 


(6) 


where  the  vector  Xt  is  defined  for  continuous  time.  The 
estimated  received  signal  correlation  matrix  is  given  by  the 
following  equation: 


nest 


i?™=- 


nest 


£  **• 


X  «  = 


nest 


■x-x1 


(7) 
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The  vector,  x^,  is  the  expression  defined  in  Equation  (6),  for 
each  discrete  time  interval.  The  quantity  nest  is  the  number 
of  estimates  used  to  compute  the  autocorrelation  matrix,  R„. 

The  autocorrelation  matrix,  R„,  may  be  used  to  extract  the 
signal  information  by  using  high  resolution  techniques  such 
as:  Multiple  Signal  Classification  (MUSIC)  [Ref.  2]  or 
Minimum  Norm  [Ref.  3] .  These  high- resolution  techniques  use 
the  singular  vector  decomposition  (SVD)  or  the  eigenvalue 
decomposition  (EVD)  of  the  received  signal  autocorrelation 
matrix  to  estimate  the  DOA  information.  A  new  SVD  (or  EVD) 
decomposition  is  needed  for  each  update  when  tracking  moving 
sources.  Since  the  SVD  (or  EVD)  decomposition  is 
computationally  expensive,  good  results  are  not  obtainable  for 
real  time  systems.  The  goal  of  this  work  is  to  find  a  good 
approximation  of  the  signal  and  the  noise  subspaces  and  to  be 
able  to  use  these  approximations  to  find  the  desired  signal 
information.  To  allow  their  use  in  real-time  algorithms,  it 
is  anticipated  that  these  approximations  will  be  a  compromise 
between  accuracy  and  processing  speed. 

The  method  to  be  used  is  the  Rank  Revealing  QR 
Factorization  (RRQR)  .  Assuming  that  AII=QR  is  the  classical  QR 
factorization  of  A,  where  Q  is  orthonormal,  R  is  upper 
triangular,  and  II  is  an  orthonormal  permutation  matrix,  this 
method  provides  a  matrix  of  the  form 
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where  the  norm  two  of  R22,  || R22 II 2 #  small  when  A  is  ill- 
conditioned.  This  is  accomplished  by  means  of  a  special 
pivoting  scheme  following  the  cited  QR  factorization.  From 
the  dimension  of  the  block  R22  of  this  matrix,  its  rank  may  be 
deduced.  This  method  also  provides  the  means  to  find 
approximations  for  the  signal  and  noise  subspaces  of  the 
matrix  A.  It  will  be  shown  that  this  information  is  contained 
in  the  orthogonal  matrix  Q,  generated  from  the  RRQR.  Details 
about  the  method  follow  in  Section  A  of  Chapter  2. 

The  RRQR  will  be  used  to  decompose  the  noise- free 
autocorrelation  matrix.  This  matrix  is  obtained  by 
subtracting  a2I  from  the  signal  autocorrelation  matrix  defined 
in  (7)  ,  where  a2  is  the  Gaussian  white  noise  variance  and  I  is 
the  identity  matrix.  The  need  for  the  noise- free 
autocorrelation  matrix  is  due  to  the  special  pivoting  scheme 
developed  by  Chan  [Ref.  4],  This  scheme  works  for  rank 
deficient  or  ill  conditioned  matrices.  However,  the 
autocorrelation  matrix,  R„,  does  not  present  this 
characteristic.  Therefore,  no  accurate  signal  subspace 
information  can  be  obtained  from  the  RRQR  decomposition  of  R„. 
To  correct  this  deficiency  of  the  algorithm,  information  about 
the  Signal  to  Noise  Ratio  (SNR)  is  needed. 
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It  is  possible  to  update  the  classical  QR  decomposition  of 
the  noise- free  autocorrelation  matrix.  The  possibility  of 
updating  the  RRQR  will  also  be  investigated,  as  this  technique 
allows  tracking  of  moving  targets.  The  rank-one  modification 
of  a  matrix  will  be  used  to  update  the  matrices  Q  and  R 
without  accessing  the  updated  noise- free  autocorrelation 
matrix. 
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II.  THE  RANK -REVEALING  QR  (RRQR)  FACTORIZATION  ALGORITHM 

This  chapter  introduces  the  theoretical  background 
necessary  to  understand  the  RRQR  algorithm  presented  by  Chan 
[Ref.  4] .  This  decomposition  will  then  be  used  to  estimate 
the  signal  information. 

A.  THE  RANK -REVEALING  QR  ALGORITHM 

The  QR  factorization  is  a  decomposition  of  a  given  m  by  n 
matrix  A  into  three  other  matrices,  Q,  R  and  II,  so  that  AlI=QR. 
The  matrix,  II  e  R““,  is  an  orthonormal  permutation  matrix,  Q 
e  Cmj“  is  orthonormal,  i.e.,  QHQ=I„,  and  R  e  C““  is  upper 
triangular.  The  QR  decomposition,  using  a  fixed  permutation 
II,  provides  a  unique  pair  of  matrices  Q  and  R.  [Ref.  4] 

When  A  is  full  rank,  R  is  non- singular .  However,  if  A  is 
rank  deficient  with  rank  deficiency  r,  we  may  choose  II  so  that 
the  rank  deficiency  of  R  is  exhibited  in  the  form  of  a  small 
lower  right  block  R22*0,  as  shown  in  Equation  (9)  .  Note  that 
since  Q  and  II  are  orthonormal  matrices,  their  singular  values 
are  equal  to  one.  Therefore,  A  and  R  have  the  same  rank. 
The  matrix  R  may  be  defined  as 

(9) 

where  R22  is  r  by  r  [Ref.  4]  . 


R= 


^11  ^121 


0  R. 


'221 
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Let  us  assume  that  <J;  is  the  i*  singular  value  of  A,  where 
<r,acr22. .  .  2<rn.  It  is  possible  to  show  that  £T0_r+1  (A)  s||r22||2  [Ref. 
4]  .  Thus,  if  the  norm  two  of  the  matrix  R22,  || R22 II 2 »  is  small, 
A  has  at  least  r  small  singular  values.  The  converse  is  not 
true,  that  is,  it  is  not  guaranteed  that,  if  A  has  r  small 
singular  values,  its  QR  factorization  will  yield  a  small 

II  ^22 II 2  • 

Thus,  it  is  necessary  to  find  an  algorithm  that  is  able  to 
reveal  the  matrix  rank  via  a  small  R22  block.  The  Rank- 
Revealing  QR  (RRQR)  factorization  algorithm  [Ref.  4]  provides 
such  a  decomposition  of  A. 

In  this  section,  a  brief  discussion  and  the  theoretical 
background  of  the  method  is  provided.  The  main  purpose  here 
is  to  identify  the  characteristics  which  are  important  to  the 
Direction  of  Arrival  (DOA)  problem.  A  thorough  theoretical 
treatment  of  the  RRQR  factorization  has  been  written  by  Chan 
[Ref.  4] . 

This  proof  starts  from  a  rank- one  deficient  matrix.  It 
will  be  extended  later  for  a  rank-r  deficient  matrix.  Let  us 
assume  that  a  given  matrix  A  is  nearly  rank-one  deficient. 
The  Theorem  2.1  of  [Ref.  4]  is  reproduced  here,  adapted  to  the 
complex  case. 

THEOREM  2 . 1  Suppose  that  we  have  a  vector  x  e  C°  with 
||x [| 2=1  such  that  ]|Ax||2=e,  and  let  II  be  a  permutation  matrix 
such  that  if  lFx=y,  then  |y„|  =  ||y||w.  Then  if  An=QR  is  the 
QR  factorization  of  All,  then 

|  rm  |  s  Vn  e  . 


8 


Proof:  First  we  note  that  since  |yn|= 
we  have 


and  ||y||2=  ||x||2=l. 


(10) 


Next,  we  have 


Q  H Ax  =  QhA IinTx  =  Ry  = 


rnnY. r\ 


(id 


Therefore, 

e  =  ||Ak|2  =  \\QhAx\\2  =  ||Ry||2>  (12) 

and  from  (10)  and  (12)  we  have  the  desired  relationship.* 

Now,  assume  that  vn  e  C°  with  ||vn||2=l  is  the  smallest  right 
singular  vector  of  A,  then  we  have 

I|AV||2=0„.  (13) 


If  we  define 

|(IPV)  |n=IML,  (14) 

where  II  is  the  permutation  as  defined  in  Theorem  2.1,  we  see 
that  All  has  a  QR  factorization  with  pivot  r,,,,  and  that  the 
magnitude  of  rm  is  less  or  equal  to  |  |  .  In  other  words, 

the  (n,n)*  element  of  R  is  small.  Therefore,  it  is  possible 
to  make  r,,,,  small  with  an  adequate  choice  of  the  permutation 
matrix. 
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Since  we  need  only  an  adequate  permutation  matrix  II,  an 
approximation  of  the  SVD  of  A  for  the  smallest  singular 
vector,  v,  is  adequate.  This  point  shows  the  advantage  of  the 
RRQR  algorithm  over  the  SVD.  An  approximation  of  v  may  be 
computed  much  faster  than  the  true  singular  vector.  In 
Chapter  3,  methods  to  find  a  good  approximation  for  v  is 
demonstrated. 

The  above  results  serve  as  the  foundation  to  compute  any 
QR  factorization  of  A.  An  approximation  of  the  smallest  right 
singular  vector  of  A  is  found.  Then  we  determine  II,  as  in 
(14)  ,  and  compute  the  QR  factorization  of  All.  Alternatively, 
the  eigenvector  corresponding  to  the  eigenvalue  closer  to  zero 
[Ref.  5]  may  be  found,  rather  than  the  smallest  right  singular 
vector  v.  The  proof  that  this  algorithm  is  valid  for  the 
eigenvector  case  is  trivial,  as  a  is  replaced  by  |\| 
throughout  the  proof. 

The  above  one -dimensional  algorithm  may  be  extended  to  the 
case  when  A  is  nearly  rank-r  deficient,  with  r>l.  We  want  to 
find  a  permutation  II,  such  that 

*u  *■ 

A"-CR*£fo  R 

is  the  QR  factorization  of  A.  The  submatrix  R^  is  (r  x  r) 
dimensional  and  J{ R22 II 2  is  small.  We  apply  the  one  dimensional 
algorithm  presented  above  to  Rn,  for  r=l,2,  ...,  where  Rn  is 


10 


the  leading  principal  (n-r  x  n-r)  dimensional  submatrix  of  R. 
After  isolating  a  (r  x  r)  dimensional  R22  block,  we  may  use 
this  one -dimensional  algorithm  to  compute  a  permutation  P  such 
that  R,,P=Q,RU  is  the  QR  factorization  of  R,,P.  Next,  we 
isolate  a  (r+1  x  r+1)  block.  This  guarantees  the  (n-r, n-r)* 
element  of  R  is  small,  leading  to: 


where 


11 

0 


Ox  %2, 


22 


and 


(16) 


(17) 


(18) 


We  notice  that  Equation  (16)  is  the  QR  factorization  of  Aft. 

The  complete  algorithm  is  summarized  below  in  steps  0  to  9. 

An  implementation  of  this  subroutine  using  the  MATLAB™ 

software  may  be  found  in  Appendix  A. 

0.  Compute  a  first  QR  decomposition  of  the  noise- free 
autocorrelation  matrix. 

"For  i=n,n-l, . . . ,n-r+l,  do: 

1.  Let  R„  be  the  leading  i  x  i  block  of  R. 

2.  Compute  the  singular  vector  v  e  C'  of  Ru,  corresponding 
to  the  minimum  singular  value  amin(Rll)  with  ||v||2=l. 
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3.  Compute  a  permutation  P  e  C1"  such  that  |  ( P“v )  1 1  =  || pnv |j . 
(This  means  find  the  maximum  absolute  value  element  of  the 
smallest  right  singular  vector  and  swap  it  with  the  i* 
element  of  the  same  vector) . 


4.  Assign  v» 


v 

0 


€  Cn  to  the  i*  column  of  W. 


5. 


Let  W=PHW,  where  P» 


P  0 
0  I ' 


6.  Compute  Q,R„,  the  QR  factorization  of  R,,P. 

7.  Let  IMIP. 


8. 


Let  Q=Q 


9.  Let  R= 


OiTP12 
0  R22 


For  this  algorithm  to  provide  the  desired  RM  block  of  R 
small  in  norm,  we  must  have  a  Rn  block  at  step  two  with  a 
small  singular  value  to  insure  that  the  (i,i)*  element  of  Ru 
is  small.  At  step  nine,  the  (n-i)*  (last)  row  of  QHiRn  must 
be  small  in  norm.  If  these  two  assertions  are  true,  then  the 
lower  (n-i+1  x  n-i+1)  block  K-22  in  step  9  is  small  in  norm  and 
the  desired  QR  factorization  exists.  To  prove  the  first 
assertion,  we  reproduce  the  Lemma  3.1  of  [Ref.  4].  Chan's 
derivation  for  the  case  of  a  real  valued  matrix  is  adapted 
here  for  the  complex  matrix  case. 
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"Lemma  3.1:  Let  B  e  C be  a  matrix  containing  any  subset 
of  k  columns  of  A.  Then 

O’min  (B)  -  <Tk(B)  S  <rk(A) 

The  submatrix  Ru  is  the  subset  of  QHAII,  composed  of  its 
first  i  columns.  QH  and  II  are  orthonormal  and  therefore  their 
singular  values  are  one.  Using  this  observation  and  the  Lemma 
3.1  of  [Ref.  4],  we  have 

°min  (R,l )  <  aj{QHAn)=(Ti(A)  . 

The  above  inequality  guarantees  that  Ru  has  a  small  singular 
value  if  (A)  is  small. 

To  prove  that  R22  is  small  in  norm,  we  reproduce  here  the 
Lemma  3.2  and  Theorems  3.1  and  3.2  of  [Ref.  4], 


Lemma  3.2:  The  matrix  W  =[wn.r+,,  .  .  . ,  wj  e  C"xr  computed  by 
the  algorithm  RRQR  satisfies  the  following  properties:  For 
i=n, n-1, . . . ,n-r+l, 


1)  |w,h-l, 

2)  (w^ j=0  for  j >i, 

3)  |  (Wi)il- 1| Wj!.  a  1/VT, 

4)  lAlfyl^  s  ci  (A)  ,  where  (Ru)  )s. 


Theorem  3.1:  Let  the  matrix  W  e  C™  as  computed  by  the 
algorithm  RRQR  be  partitioned  as  W“=  (W,H, W2H)  ,  where  W*1  e 
C™  is  upper  triangular  and  non- singular .  Then  the  QR 
factorization  of  All  as  computed  by  the  algorithm  RRQR 
satisfies : 


R22I2  s  ^n-r+l  W2  ‘  1 2  Vr. 
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Proof:  Denote  the  columns  of  W  by  {wn.r+l,  .  .  . ,  wn} .  Define 
Y*QHAlIW  €  C"”.  From  property  (4)  of  Lemma  3.2,  we  have 


llil 


2  ^  114  ~ 


n-x*l 


*  vT7T6n.r+1  s  /TiTon  -r+1 ' 


(23) 


where  6;  »  (a^ (Ru)  )  {  and  ||.|jF  denotes  the  Frobenius  norm.  Next 
the  matrix  Y  can  be  expressed  as : 


Y  =  Q  tARW  =  RW= 


Rl2^2 


(24) 


which  leads  to 


||  Y\\2-z\\R22W2\\2± 


\R22  »2 
w2-%  ‘ 


(25) 


Combining  Equations  (23) ,  (24)  and  (25) ,  we  get 


Il*22l2 
IIWz'1!  zS*' 


(26) 


from  which  the  desired  result  follows.* 

Theorem  3.2:  The  algorithm  RRQR  computes  a  permutation 
II  and  a  QR  factorization  of  A,  given  by  AII=QR  where  the 
elements  of  the  lower  (r  x  r)  upper  triangular  block  of  R 
satisfy 

\rU I  *  okyfl c  . 

k-i  •  ' 

s  2j~ior/n  for  n-r<izj  zn 


Proof:  Using  Lemma  3.2,  we  have,  for  n-rcisj  s  n 
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a^\AIlw^2=lRWji2i\  J  rlk(wj)  k\ . 

k-i 


Isolating  the  k=j  term  in  the  sum  and  rearranging  terms,  we 
get 

j-i 

I  *  °j+E  *1  • 

k=i 


From  Lemma  3.2,  |  (Wj)  i  \  =  ||  (Wj)  H.  a  I/n/J,  we  have 

kijl  * 

k~i 


Solving  this  recurrence  in  the  index  j ,  we  get  the  bound  given 
in  the  first  inequality  in  Equation  27.  Using  the  bound  o$A c 
s  cry n  in  each  term  of  the  sum  in  the  desired  result,  we  get 
the  second  bound.* 

Equation  (27)  in  Theorem  3.2  shows  the  bounds  for  the 
elements  of  the  matrix  R  generated  by  the  RRQR  algorithm.  The 
second  inequality  states  the  bounds  for  the  elements  of  the 
block  R22.  The  factor  2y'  indicates  that  one  element  in  RM 
increases  with  its  distance  from  the  main  diagonal.  There¬ 
fore,  for  large  values  of  the  rank  deficiency  r,  the  bounds 
may  be  quite  large,  as  it  grows  exponentially.  Numerical 
simulation  cases  will  be  shown  later  to  demonstrate  that  these 
bounds  are  overconservative. 

This  algorithm  was  originally  proposed  to  be  iterated  from 
n  until  n-r+1,  where  n  is  the  number  of  sensors  in  the  array 
and  r  is  the  noise- free  autocorrelation  matrix  rank  deficien¬ 
cy.  Alternatively,  Prasad  and  Chandna  [Ref.  5]  proposed  to 
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iterate  the  recursion  until  n-m+1,  where  m  is  the  number  of 
sources.  A  faster  algorithm  is  the  result  of  this  technique 
as  generally  there  are  fewer  sources  than  sensors.  The 
drawback  of  this  approach,  however,  is  a  loss  of  precision  in 
the  estimated  signal  information,  as  the  algorithm  might  not 
capture  all  the  rank  deficiency  of  the  matrix  to  be  decom¬ 
posed.  This  problem  will  be  investigated  in  Chapter  4. 

B.  USING  THE  RRQR  ALGORITHM  TO  FIND  THE  NOISE  AND  SIGNAL 

SUBSPACES  OF  THE  NOISE -FREE  AUTOCORRELATION  MATRIX 

Assume  that  we  have  a  linear  equispaced  array  composed  of 
n  sensors  and  m  sources,  with  n  >  m.  The  theoretical  noise 
free  autocorrelation  matrix,  (R.)  ,  is  r  (=n-m)  rank  deficient. 
However,  the  practical  noise- free  autocorrelation  matrix  is 
nearly  rank  deficient. 

It  is  possible  to  find  the  signal  and  noise  subspaces  of 
an  incoming  signal  using  the  RRQR  algorithm  shown  in  Section 
A.  After  applying  the  complete  RRQR  algorithm,  the  matrix  R 
reveals  the  rank  deficiency  of  R,.  Consequently,  the  norm  of 
R22  is  small  compared  to  the  norm  of  the  rest  of  the  matrix  R. 
Note  that  R  is  upper  triangular.  Here,  Rn  has  dimension  (n- 
r(=m)  x  n-r)  ,  R12  has  dimension  ((n-r)  x  r)  ,  and  R^  has 
dimension  (r  x  r) .  The  signal  subspace  is  contained  in  the 
first  m  columns,  and  the  noise  subspace  is  contained  in  the  r 
last  columns  of  the  matrix  Q.  The  matrix,  W,  of  the  RRQR 
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algorithm  is  not  used  here  as  an  approximation  of  the  noise 
subspace,  as  originally  proposed  by  Chan  [Ref.  4]  .  The  use  of 
this  matrix  to  approximate  the  noise  subspace  yielded  poor 
results,  as  shown  by  Fargues  [Ref.  3]. 


C.  USING  THE  GIVENS  ROTATION  TO  REPACTOR  R„P 

A  new  QR  factorization  of  the  RnP  term  is  presented  in 
step  six  of  the  RRQR  algorithm.  In  step  one,  it  is  reasonable 
to  run  a  complete  QR  algorithm  for  the  first  QR  factorization. 
But  as  Ru  is  already  upper  triangular  and  P  is  only  a  permuta¬ 
tion  that  changes  the  position  of  two  rows,  we  can  use 
inclusive  Givens  Rotations  to  zero  out  the  few  non- zero 
elements  of  RnP  below  its  main  diagonal .  These  rotations 
yield  a  more  efficient  algorithm  to  deal  with  the  special 
structure  of  the  problem. 

The  complex  Givens  Rotation  matrix  is  defined  as 

c  s 

G=  X 

-conj(s)  c 


(31) 


where  c  is  real,  s  is  complex  and  such  that  |c|2+|s|2=l. 

The  Givens  Rotation  matrix  is  defined  such  that  the 
following  equality  is  satisfied  [Ref.  6] : 


G* 


X 

i 

y 

0 

(32) 


Suppose  we  have  a  matrix,  A,  and  we  want  to  zero  out  a 
given  element  a(k,i),  below  the  main  diagonal,  i.e.,  k>i.  If 
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we  multiply  this  matrix  by  the  matrix  J  defined  below,  we  will 
make  the  element  a(k,i)  equal  to  zero.  If  we  repeat  this 
procedure  for  all  elements  of  A  below  the  main  diagonal  and 
different  from  zero,  the  matrix  A  is  transformed  into  an  upper 
diagonal  matrix.  This  is  the  essence  of  the  Givens  Rotation. 


v  o  o 

o  c  o  $  o 

:  "v 

0  0  0  0 

1 

0  -<x>n)(t>  0  ft  0 

0  0  1 


I 


k 


i  k 

The  subroutine  GIVENS1  used  to  compute  the  complex  Givens 
Rotation  matrix  J  and  the  subroutine  QRGIV  used  to  perform  the 
QR  factorization  using  the  Givens  rotation  may  be  found  in 
Appendices  B  and  C. 
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III.  METHODOLOGY  FOR  APPROXIMATING  THE  SMALLEST  RIGHT 
SINGULAR  VECTOR 

Chapter  2  presented  the  RRQR  algorithm  and  indicated  that 
a  procedure  is  needed  to  find  reliable  approximations  of  the 
smallest  right  singular  vector  of  a  matrix.  The  inverse 
iteration  method  and  the  Incremental  Condition  Estimator  are 
introduced  to  address  this  problem.  The  inverse  iteration 
method  is  used  to  find  the  smallest  singular  vectors  of  a 
matrix.  It  may  also  be  used  to  find  the  smallest  eigenvector 
in  absolute  value.  The  Incremental  Condition  Estimator  is 
used  to  find  a  reliable  approximation  of  the  smallest  singular 
value  and  the  corresponding  singular  vectors  of  a  matrix. 

A.  THE  INVERSE  ITERATION  METHOD  TO  FIND  THE  SMALLEST  SINGULAR 

VECTOR 

Step  two  of  the  RRQR  algorithm  is  a  computationally  expen¬ 
sive  part  of  the  procedure  [Ref.  7] .  It  is  therefore  desir¬ 
able  to  find  an  algorithm  that  finds  the  minimum  singular 
value  and  its  corresponding  right  singular  vector  quickly  and 
accurately. 

The  algorithm  implemented  is  the  inverse  iteration  method. 

Starting  with  an  initial  guess,  v0  (in  this  case  a  vector 
composed  of  ones) ,  the  following  algorithm  is  iterated  until 
convergence  [Ref.  8] . 
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1.  Solve  Aui+l=Vj  for  ui+1, 

2.  Let  ui+1=ui+1/||ui+1||2, 

3.  Solve  AH^i+1-ui+,  for  vj+l, 

4.  Let  vi+1=vi+,/||vi+,||2. 

Let  the  converged  be  v,,,  and  compute  and  by: 

Au=vw,  u„»u/||u||2,  amin=l/ ||a||2. 

The  above  algorithm  computes  the  right  (v„)  and  left  (u„) 
singular  vectors  as  well  the  minimum  singular  value  of  the 
matrix  A.  This  algorithm  yields  good  results  for  a  small 
number  of  iterations,  as  shown  by  the  numerical  simulations 
presented  below.  Three  iterations  were  used  here  to  estimate 
the  singular  vectors. 

To  evaluate  the  results  achieved  with  the  inverse  itera¬ 
tion  algorithm,  we  examined  four  different  test  cases.  One  to 
three  iterations  were  run  on  each  test  case.  In  each  case, 
the  right  singular  vector  was  compared  to  the  corresponding 
vector  generated  by  the  SVD  decomposition  computed  with  the 
MATLAB™  software.  This  procedure  was  followed  for  1000 
different  random  matrices. 

Comparisons  between  approximated  and  computed  singular 
vectors  were  obtained  by  evaluating  the  magnitude  of  the 
projection  of  the  estimate  over  the  true  smallest  right 
singular  vector.  If  the  two  vectors  are  parallel,  the 
projection  is  maximum  and  equal  to  one.  If  they  are  perpen¬ 
dicular,  the  projection  is  zero.  The  projections  were 
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distributed  among  ten  bins,  ranging  from  0  to  1,  as  shown  in 
Tables  1-4.  The  mean  and  standard  deviation  for  each  one  of 
the  iteration  steps  were  computed. 

The  first  test  case  used  a  10x10  dimensional  random  matrix 
generated  using  MATLAB™.  The  second  one  used  a  10x10  random 
matrix,  with  the  lower  triangle  was  imposed  as  zero.  The 
third  test  case,  used  an  upper  triangular  matrix  R,  obtained 
from  the  QR  decomposition  without  pivoting  of  a  10x10  random 
matrix.  Finally,  the  fourth  test  case  used  an  upper  triangu¬ 
lar  matrix  R  obtained  from  the  QR  decomposition  with  pivoting 
of  a  10x10  random  matrix.  These  cases  were  evaluated  in  terms 
of  performance  to  verify  the  adequacy  of  the  method  within 
different  and  perhaps  more  demanding  contexts.  The  test 
results  are  summarized  in  Tables  1  through  5. 

Table  1:  INNER  PRODUCT  MAGNITUDE  BETWEEN  ESTIMATED  AND  TRUE 
SMALLEST  SING.  VECTORS,  RANDOM  MATRIX,  1000  TRIALS. 


Percentage 

of  Projec¬ 
tion  lying 
between 

One 

Itera 

tion 

Three 

Itera 

tions 

0  and 

.i 

.6 

.i 

.1 

.1  and 

.2 

1.4 

.6 

.2 

.2  and 

.3 

.7 

.2 

.1 

.3  and 

.4 

.5 

.7 

.4 

.4  and 

.5 

.8 

.7 

.1 

.5  and 

.  6 

.9 

.4 

.3 

.  6  and 

.7 

1.1 

0 

.4 

.7  and 

.8 

2.4 

1.0 

.6 

.  8  and 

.9 

3.9 

.7 

.8 

.  9  and 

1 

87.7 

95.6 

97.0 

.9441 

.9756 

.9867 

5535S3I 

.1645 

.1121 

.0802 
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Table  2:  INNER  PRODUCT  MAGNITUDE  BETWEEN  ESTIMATED  AND  TRUE 
SMALLEST  SING.  VECTOR,  RANDOM  MATRIX  WITH  LOWER  TRIANGLE 
PORTION  EQUAL  TO  ZERO,  1000  TRIALS. 


Percentage 
of  Projec¬ 
tion  lying 
between 

One 

Itera 

tion 

Two 

Itera 

tions 

Three 

Itera 

tions 

0  and  . 1 

0 

0 

0 

. 1  and  . 2 

0 

0 

0 

.2  and  .3 

.1 

0 

0 

.3  and  .4 

.1 

0 

0 

. 4  and  . 5 

0 

.1 

0 

. 5  and  . 6 

.1 

0 

0 

. 6  and  . 7 

.  1 

0 

0 

. 7  and  . 8 

.3 

0 

0 

. 8  and  . 9 

.5 

.4 

.1 

.9  and  1 

98.8 

99.5 

99.9 

Average 

.9960 

.9988 

.9996 

Standard 

Deviation 

.0399 

.0184 

.0071 

Table  3:  INNER  PRODUCT  MAGNITUDE  BETWEEN  ESTIMATED  AND  TRUE 
SMALLEST  SING.  VECTORS,  UPPER  TRIANGULAR  MATRIX  RESULTING  FROM 
QR  FACT.  W/O  PIVOTING,  1000  TRIALS. 


Percentage  of 
Projection 
lying  between 

One 

Iterat 

ion 

Two 

Iterat 

ions 

Three 

Iterati 

ons 

0  and  . 1 

.7 

.2 

.2 

. 1  and  . 2 

.5 

.6 

.1 

. 2  and  . 3 

.7 

.2 

.4 

.3  and  .4 

.  6 

.5 

.4 

. 4  and  . 5 

1.2 

0 

.1 

. 5  and  . 6 

.6 

.1 

.1 

. 6  and  . 7 

.9 

.7 

.1 

. 7  and  . 8 

2.0 

1.1 

.2 

.8  and  .9 

2.8 

1.4 

.9 

.9  and  1 

90.0 

95.2 

97.5 

.9528 

.9781 

.9868 

mum 

.1480 

.1044 

.0843 
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Table  4:  INNER  PRODUCT  MAGNITUDE  BETWEEN  ESTIMATED  AND  TRUE 
SMALLEST  SING.  VECTORS,  UPPER  TRIANGULAR  MATRIX  RESULTING  FROM 
QR  FACT.  W/  PIVOTING,  1000  TRIALS. 


Percentage 
of  Projec¬ 
tion  lying 
between 

One 

Itera 

tion 

TWO 

Itera 

tions 

Three 

Itera 

tions 

0  and  . 1 

0.1 

.1 

.1 

.1  and  .2 

.6 

.3 

0 

. 2  and  . 3 

.5 

.1 

.3 

. 3  and  . 4 

.7 

.2 

.1 

. 4  and  . 5 

.7 

.  6 

.2 

. 5  and  . 6 

.7 

.6 

.4 

. 6  and  . 7 

.4 

.6 

.5 

. 7  and  . 8 

.7 

.5 

.4 

. 8  and  . 9 

3.0 

.1 

.6 

.9  and  1 

92.6 

96.9 

97.4 

Average 

.9660 

.9823 

.9887 

Standard 

Deviation 

.1225 

.0922 

.0743 

Table  5:  COMPARISON  BETWEEN  THE  FOUR  CASES  EXAMINED  IN  TABLES 
1  THROUGH  4. 


Totally 

Random 

Matrix 

Imposed 

upper 

triangul 

ar 

Resulting 

from  QR 
decomp, 
w/o  pivot - 
ing 

Resulting 

from  QR 
decomp .  w/ 
pivoting 

Percentage  of 
projections 
between  . 9  and 
1.0 

97 

99.9 

97.5 

97.4 

.9867 

.9996 

.9868 

.9887 

Standard  Devi¬ 
ation  of  Pro¬ 
jections 

.0802 

.0071 

.0843 

.0743 
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Table  5  shows  that  the  best  result  was  obtained  from  the 


imposed  triangular  matrix  experiment.  However,  all  of  the 
experiments  show  good  results.  These  results  justify  the  use 
of  the  inverse  iteration  method  as  a  low  cost  alternative  to 
approximate  the  smallest  right  singular  vector  in  step  two  of 
the  RRQR  algorithm.  The  important  point  is  to  find  a  permuta¬ 
tion  so  that  the  element  of  the  least  dominant  singular  vector 
that  presents  the  largest  magnitude  be  positioned  at  the  i* 
row  [Ref.  4]  .  Therefore,  it  is  not  necessary  to  find  the 
exact  singular  vector  for  the  RRQR  algorithm.  The  source  code 
implemented  to  approximate  the  smallest  singular  vector  by  the 
use  of  inverse  iteration  method  is  shown  in  Appendix  J. 

B.  THE  INVERSE  ITERATION  METHOD  TO  FIND  THE  SMALLEST 

EIGENVECTOR 

Prasad  [Ref.  5]  states  that  the  RRQR-based  algorithm  works 
when  using  an  EVD  (Eigenvector  Decomposition)  rather  than  a 
SVD  (Singular  Vector  Decomposition)  of  the  noise- free 
autocorrelation  matrix.  The  validity  of  this  approach  is 
investigated  here.  Theoretically,  it  is  correct  to  only 
replace  ff  by  |A|  in  the  RRQR  factorization  proof  presented  in 
Chapter  2 . 

The  inverse  iteration  method  may  be  used  to  find  the  least 
dominant  eigenvector  (the  one  with  a  magnitude  closer  to 
zero)  .  Again,  we  find  an  initial  guess  v0,  which  may  be  a 
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vector  composed  solely  of  ones.  The  following  algorithm 
should  be  iterated  until  convergence  of  v. 

1.  Solve  Ayk=vk.,,  for  yk. 

2.  Normalize  vk=yk/ ||yk||2. 

Comparing  the  two  algorithms,  we  note  that  this  inverse 
iteration  algorithm  uses  the  same  initial  two  steps  as  the  one 
used  for  the  minimum  singular  value.  To  compare  the 
performance  of  these  two  algorithms,  the  same  test  cases  were 
run  as  in  previous  section.  However,  six  iterations  were  run 
to  provide  a  viable  comparison.  The  results  are  shown  in 
Tables  6  through  9 . 
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Table  6:  INNER  PRODUCT  MAGNITUDE  BETWEEN  ESTIMATED  AND  TRUE 
SMALLEST  EIGENVECTORS ,  RANDOM  MATRIX,  1000  TRIALS. 


Percentage 

of  Projec¬ 
tion  lying 
between 

H 

Three 

Iterati 

ons 

Four 

Iterati 

ons 

Five 

Iterati 

ons 

Six 

Iterat 

ions 

0  and  . 1 

4.0 

2.1 

2.0 

0.6 

an 

. 1  and  . 2 

4.4 

2.3 

1.3 

1.4 

1.3 

1.2 

. 2  and  . 3 

3.6 

2.6 

1.7 

1.4 

1.8 

1.0  ! 

. 3  and  . 4 

5.3 

3.4 

2.5 

2.3 

1.4 

. 4  and  . 5 

6.8 

4.0 

EBm 

2.6 

2.8 

an 

. 5  and  . 6 

an 

6.7 

5.7 

4.1 

3.4 

mam 

. 6  and  . 7 

9.3 

7.0 

5.0 

6.8 

4.7 

^0 _ II 

. 7  and  . 8 

11.2 

11.1 

7.5 

6.7 

7.7 

ran 

. 8  and  . 9 

17.2 

13.0 

10.6 

11.0 

9.6 

9.3 

.9  and  1 

30.9 

47.8 

59.7 

63.1 

66.4 

66.7 

HW  A  t»  /-—MB 

.6964 

.7861 

.8353 

.8586 

.8694 

.8775 

.2758 

.2491 

.2347 

.2124 

.2116 

.1986 

Table  7  :  INNER  PRODUCT  MAGNITUDE  BETWEEN  ESTIMATED  AND  TRUE 
SMALLEST  EIGENVECTORS,  RANDOM  MATRIX  WITH  LOWER  TRIANGLE 
PORTION  EQUAL  TO  ZERO,  1000  TRIALS. 


Percentage 
of  Projec¬ 
tion  lying 
between 

One 

Itera 

tion 

Two 

Itera 

tions 

Three 

Itera 

tions 

Four 

Itera 

tions 

Five 

Itera 

tions 

Six 

Itera 

tions 

0  and  . 1 

0.9 

0.5 

0.4 

0.3 

0.5 

0.3 

. 1  and  . 2 

0.7 

0.7 

0.5 

0.3 

0.4 

0.2 

. 2  and  . 3 

1.5 

0.6 

0.4 

0.2 

0.1 

0.2 

.3  and  .4 

0.7 

0.6 

0.4 

0.4 

0.4 

0.2 

.4  and  .5 

1.3 

1.0 

0.7 

1.1 

0.5 

0.7 

. 5  and  . 6 

1.9 

1.5 

0.6 

0.9 

0.5 

0.3 

. 6  and  . 7 

1.4 

1.4 

1.1 

1.1 

0.4 

0.9 

. 7  and  . 8 

3.5 

2.9 

1.1 

0.8 

0.9 

0.4 

. 8  and  . 9 

5.7 

4.2 

2.5 

2.6 

1.7 

2.1 

.9  and  1 

82.4 

86.6 

92.3 

92.3 

94.6 

94.7 

Average 

.9217 

.9445 

.9634 

.9662 

.9740 

.9774 

Standard 

Deviation 

.1787 

.1473 

.1280 

.1201 

.1134 

.1019 
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Table  6:  INNER  PRODUCT  MAGNITUDE  BETWEEN  ESTIMATED  AND  TRUE 
SMALLEST  EIGENVECTOR,  UPPER  TRIANGULAR  MATRIX  RESULTING  PROM 
QR  PACT.  W/O  PIVOTING,  1000  TRIALS. 


Percentage  of 
Projection 
lying  between 

One 

Iterat 

ion 

Three 

Iterat 

ions 

Four 

Iterat 

ions 

Five 

Iterat 

ions 

Six 

Iterat 

ions 

0  and  . 1 

3.1 

2.0 

0.8 

1.0 

0.3 

1.2 

. 1  and  . 2 

2.3 

1.0 

1.5 

1.4 

1.4 

0.9 

. 2  and  . 3 

3.1 

1.3 

1.0 

1.2 

0.9 

0.7 

. 3  and  . 4 

3.3 

1.8 

1.5 

1.2 

1.3 

0.6 

. 4  and  . 5 

3.7 

1.8 

1.9 

1.3 

2.0 

. 5  and  . 6 

6.1 

3.5 

1.7 

1.8 

1.2 

1.0 

. 6  and  . 7 

6.9 

4.1 

3.1 

1.1 

1.5 

0.8 

. 7  and  . 8 

8.9 

6.6 

4.4 

2.5 

2.3 

1.8 

. 8  and  . 9 

17.3 

10.6 

6.0 

6.0 

3.9 

3.2 

.9  and  1 

44.1 

65.4 

78.2 

81.9 

85.9 

00 

• 

r*- 

CO 

.7654 

.8562 

.9018 

.9158 

.9328 

.9379 

.2579 

.2227 

.1943 

.1945 

.1728 

.1767 

Table  9:  INNER  PRODUCT  MAGNITUDE  BETWEEN  ESTIMATED  AND  TRUE 
SMALLEST  EIGENVECTORS,  UPPER  TRIANGULAR  MATRIX  RESULTING  FROM 
QR  FACT.  W/  PIVOTING,  1000  TRIALS. 


Percentage 
of  Projec¬ 
tion  lying 
between 

One 

Itera 

tion 

Two 

Itera 

tions 

Three 

Itera 

tions 

Four 

Itera 

tions 

Five 

Itera 

tions 

Six 

Itera 

tions 

0  and  . 1 

1.3 

0.4 

0.3 

0.0 

0.2 

0.1 

. l  and  . 2 

0.8 

0.7 

0.4 

0.4 

0.0 

0.1 

. 2  and  . 3 

1.1 

0.5 

0.2 

0.3 

0.0 

0.0 

. 3  and  . 4 

1.5 

0.4 

0.5 

0.2 

0.4 

0.0 

.4  and  .5 

3.1 

0.8 

0.4 

0.2 

0.3 

0.2 

. 5  and  . 6 

3.7 

1.2 

0.8 

0.4 

0.3 

0.3 

. 6  and  . 7 

5.6 

2.2 

0.8 

0.7 

0.5 

0.4 

. 7  and  . 8 

9.8 

3.2 

1.7 

2.0 

1.1 

1.0 

.8  and  .9 

16.9 

7.1 

5.0 

3.1 

2.6 

2.4 

.9  and  1 

56.2 

83.5 

89.9 

92.7 

94.6 

95.5 

Average 

.8451 

.9344 

.9611 

.9734 

.9807 

.9856 

Standard 

Deviation 

.1969 

.1432 

.1167 

.0949 

.0803 

.0668 
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The  results  found  in  Tables  1  through  4  and  6  through  9 
are  summarized  in  Table  10  below.  Again,  the  case  for  imposed 
triangular  matrices  showed  the  best  performance. 


Table  10:  COMPARISON  BETWEEN  THE  FOUR  CASES  CONSIDERED  FOR  THE 
SVD  VERSUS  THE  EIGENVECTOR  INVERSE  ITERATION  METHOD. 


Totally 

Random 

Matrix 

Imposed 
upper  tri¬ 
angular 

Resulting 

from  QR 
decomp . 
w/o  pivot¬ 
ing 

Resulting 

from  QR 
decomp,  w/ 
pivoting 

SVD 

EIGV 

SVD 

EIGV 

SVD 

EIGV 

SVD 

EIGV 

Percentage  of 
projections 
between  . 9 
and  1.0 

97.0 

66.7 

99.9 

94.7 

97.5 

87.8 

97.4 

95.5 

Average  of 
projections 

.9867 

.8775 

.9996 

.9774 

.9868 

.9379 

.9887 

.9856 

Standard  De¬ 
viation  of 
Projections 

.0802 

.0071 

.1019 

.0843 

.1767 

.0743 

.0668 

The  SVD  columns  show  the  results  of  three  iterations  of 
the  inverse  iteration  method  to  find  the  singular  vector 
corresponding  to  the  minimum  singular  value.  The  columns  EIGV 
are  the  results  of  six  iterations  of  the  inverse  iteration 
method  to  find  the  least  dominant  eigenvector. 

A  cursory  look  at  these  results  show  that  the  eigenvector 
approximation  is  worse  than  those  for  the  singular  vectors 
approximation.  A  hypothesis  test  based  on  the  two  saitqples  is 
therefore  tested  for  each  case.  An  assumption  is  made  that 
each  sample  came  from  different  populations  with  each  group  of 
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1000  projections  considered  to  be  one  sample.  The  assumptions 


are 


1.  The  first  sample,  for  the  eigenvector  case,  is  a  random 
sample  from  a  population  with  mean  pt,  and  standard 
deviation  ax . 

2.  The  second  sample,  for  the  SV  case,  is  a  random  sample 
from  a  population  with  mean  fa  and  standard  deviation  a2. 

3.  Both  samples  are  independent  of  one  another. 

4.  The  samples  are  large  enough  to  apply  the  Central  Limit 
Theorem. 

The  hypothesis  test  can  be  described  by  [Ref.  9] 

Ho:  Mi=M2 

H|  2  jLij  <  fa  » 


One  hypothesis  test  will  be  conducted  for  each  one  of  the 
four  test  cases.  At  a  level  of  significance  of  1%,  Hq  will  be 
rejected  if  z  s  -2.33,  where 


avg{x)  -avg{y) 


1 


2  2 
n  n 


(33) 


If  Hq  is  false,  we  may  decide  that  fa  is  smaller  than  fa. 
In  Equation  (33),  n  is  the  sample  size  of  1000,  avg(x)  and 
avg(y)  are  the  averages  to  be  compared  and  Si  and  s2  are  the 
standard  deviations  computed  from  the  samples.  For  example, 
the  first  test  case  avg (x) = . 8775 ,  avg (y) = . 9867,  s^.1986  and 
s2=.0802  from  Table  10.  Applying  these  values  to  Equation 
(33),  we  find  z=-16.12. 
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Table  11:  RESULTS  FOR  THE  HYPOTHESIS  TEST  CONDUCTED 
ON  THE  COMPARISON  BETWEEN  THE  SVD  VERSUS  THE 
EIGENVECTOR  INVERSE  ITERATION  METHOD. 


i 

-16.12  Reject  Hq  =>  pi,  <  pi2 

2 

-6.87  Reject  Hq  =>  pi,  <  pi2 

3 

-7.90  Reject  Ho  =>  pt,  <  pi2 

4 

-.98  jCannot  Reject  Hq  j 

As  can  be  seen  in  Table  11  above,  the  first  average  pi,  is 
significantly  smaller  than  pi2  (at  1%  of  level  of  significance, 
z<-2.33),  for  all  cases  but  case  four.  Hence,  we  reject  Hq 
for  the  first  three  cases.  These  results  justify  the  choice 
of  using  the  singular  vector  approximated  by  the  inverse 
iteration  method  for  use  with  the  algorithm.  Note  that  there 
is  no  theoretical  reason  for  not  using  the  smallest  eigenvec¬ 
tor  rather  than  the  smallest  right  singular  vector  in  the  RRQR 
algorithm.  They  span  the  same  subspace.  The  reasons  for 
choosing  the  singular  vector  relies  solely  on  the  fact  that 
the  inverse  iteration  method  yields  better  results.  The 
source  code  implemented  to  approximate  the  smallest  eigenvec¬ 
tor  by  the  inverse  iteration  method  is  shown  in  Appendix  M. 

C.  THE  INCREMENTAL  CONDITION  ESTIMATOR 

A  third  method  tested  here  to  estimate  the  singular  vector 
corresponding  to  the  minimum  singular  value  is  the  Incremental 
Condition  Estimator  (ICE)  [Ref.  10]  .  Suppose  that 
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A=  [a,, . . .  ,aj  is  a  (m  x  n)  dimensional  matrix  and  let  ex,  a...  a 
Gnn  a  0  be  its  singular  values.  The  minimum  singular  value  of 
A  measures  how  close  this  matrix  is  to  rank  deficiency.  The 
condition  number  becomes 


k2(A) 

°min 


(34) 


Now,  suppose  we  take  the  QR  factorization  of  A.  This 
algorithm  is  intended  to  work  from  a  lower  triangular  matrix 
L.  Since  R  is  upper  triangular  let  us  define  L=RT.  There 
will  be  no  loss  of  generality  here  because  the  singular  values 
of  any  matrix  and  its  transpose  are  the  same  [Ref.  10] . 

If  we  have  a  n  dimensional  lower  triangular  matrix  L 
generating  one  row  at  a  time  with  an  approximate  singular 
vector  x  such  that  (L) «»1/ |jx||2  and  a  new  row  (vT,7)  of  L  ,  we 
can  obtain 


L'  * 


L  0 
VT  Y 


(35) 


such  that  amin(L'  )~l/||y||2  without  having  to  access  L  again  [Ref. 
10]  . 

Given  x  such  that  Lx=d  with  ||d||2=l  and  (L)  —1/  ||x||2,  let 
us  find  s=sin(0)  and  c=cos{$)  such  that  ||y||2  is  maximized, 
where  y  solves 


L'y  = 


L  0 
vr  Yl 


y  = 


sen 


=  d' 


(36) 
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The  solution  for  the  above  equation  is 


y  = 


sx 

c-sa  , 
Y 


where  a=vTx  and  ||d' ||2=||d||2=l . 
Now  define 


(37) 


1  +  P  -a 
a  1 


where  (3=y2xrx+a2-l .  In  this  case,  we  have 


(38) 


Ml  =  -L|s  cJB 
Y2 


s 

c 


(39) 


Assuming  y?*0 ,  the  optimal  pair  (s,c)T  is  the  eigenvector 
corresponding  to  the  largest  eigenvalue  of  B  [Ref.  10] . 
Also  assuming  as*0,  we  may  define 

and  \i=i\+sign(a)  (40) 

2tt 

to  obtain  \m„=au+l .  Finally  the  optimal  pair  (s,c)T  is  given 
as 


i  -  1 

H  Vn2+1 

-1 

(41) 


After  computing  the  optimal  pair  (s,c)T,  a  new  approximate 
singular  vector,  y,  may  be  computed  as  defined  in  Equation 
(37)  and  the  smallest  singular  value  of  L'  by 
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(42) 


®min 


(L'>« 


1 


Using  the  above  estimator  we  ran  the  same  cases  as  for  the 
inverse  iteration  presented  in  Sections  A  and  B  above.  The 
results  are  summarized  in  Tables  12  through  15. 


Table  12:  INNER  PRODUCT  MAGNITUDE  BETWEEN  ESTIMATED  AND  TRUE 
SMALLEST  SING.  VECTORS,  RANDOM  MATRIX,  1000  TRIALS. 


Percentage  of 
Proj  ection 
lying  between 

0  and  . 1 

0.6 

. 1  and  . 2 

0.2 

. 2  and  . 3 

0.5 

. 3  and  . 4 

1.1 

. 4  and  . 5 

0.7 

.5  and  .6 

0.6 

. 6  and  . 7 

0.9 

. 7  and  . 8 

. 8  and  . 9 

5.6 

.9  and  l 

87.7 

.9015 

Standard  De¬ 
viation 

.2041 
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Table  13:  INNER  PRODUCT  MAGNITUDE  BETWEEN  ESTIMATED  AND  TRUE 
SMALLEST  SING.  VECTORS,  RANDOM  MATRIX  WITH  LOWER  TRIANGLE 
PORTION  EQUAL  TO  ZERO,  1000  TRIALS. 


Percentage  of 
Projection 
lying  between 

0  and  . 1 

0.0 

. 1  and  . 2 

0.0 

. 2  and  . 3 

0.0 

. 3  and  . 4 

0.0 

.4  and  .5 

0.2 

. 5  and  . 6 

0.3 

. 6  and  . 7 

. 7  and  . 8 

0.3 

.8  and  .9 

0.5 

.9  and  1 

98.4 

Average 

.9929 

WrmXlSM 

.0441 

Table  14:  INNER  PRODUCT  MAGNITUDE  BETWEEN  ESTIMATED  AND  TRUE 
SMALLEST  SING.  VECTORS,  UPPER  TRIANGULAR  MATRIX  RESULTING  FROM 
QR  FACT.  W/O  PIVOTING,  1000  TRIALS. 


Percentage  of 
Projection 
lying  between 

0  and  . 1 

0.3 

. 1  and  . 2 

0.3 

.2  and  .3 

0.4 

.3  and  .4 

0.1 

.4  and  .5 

0.5 

. 5  and  . 6 

an 

. 6  and  . 7 

0.9 

. 7  and  . 8 

0.7 

. 8  and  . 9 

u 

. 9  and  1 

93.6 

Average 

.9705 

Standard  De¬ 
viation 

.1079 
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Table  15:  INNER  PRODUCT  MAGNITUDE  BETWEEN  ESTIMATED  AND  TRUE 
SMALLEST  SING.  VECTORS,  UPPER  TRIANGULAR  MATRIX  RESULTING  FROM 
QR  FACT.  W/  PIVOTING,  1000  TRIALS. 


Percentage  of 
Projection 
lying  between 

0  and  . 1 

0.3 

. 1  and  . 2 

0.2 

. 2  and  . 3 

0.0 

. 3  and  . 4 

0.2 

.4  and  .5 

0.1 

. 5  and  . 6 

0.1 

.6  and  .7 

0.2 

.7  and  .8 

0.7 

. 8  and  . 9 

1.3 

.9  and  1 

96.9 

Average 

.9859 

.0797 

Table  16  summarizes  the  results  found  for  each  one  of  the 
four  test  cases,  comparing  the  inverse  iteration  method  to 
find  the  smallest  right  singular  vector  and  the  Incremental 
Condition  Estimator.  As  before,  those  numbers  represent  the 
percentage  of  projections  lying  between  .9  and  1  of  the 
corresponding  vectors  found  via  the  corresponding  approxima¬ 
tion  method  over  the  true  smallest  right  singular  vector  found 
via  SVD . 


35 


Table  16:  COMPARISON  BETWEEN  THE  POUR  CASES  CONSIDERED  FOR  THE 
SVD  INVERSE  ITERATION  VERSUS  THE  INCREMENTAL  CONDITION 
ESTIMATOR  METHOD. 


Totally 

Random 

Matrix 

Imposed 
upper  tri¬ 
angular 

Resulting 
from  QR 
decomp . 
w/o  pivot¬ 
ing 

Resulting 
from  QR 
decomp,  w/ 
pivoting 

INV. 

ITER. 

ICE 

INV. 

ITER. 

ICE 

INV. 

ITER. 

ICE 

INV. 

ITER. 

ICE 

Percentage  of 
projections 
between  . 9 
and  1.0 

97.0 

87.7 

98.4 

97.5 

93.6 

97.4 

96.9 

smsfEESH 

.9867 

.9015 

.9996 

.9929 

.9868 

.9705 

.9887 

.9859 

Standard  De¬ 
viation  of 
Projections 

.0802 

.2041 

.0071 

.0441 

.0843 

.1079 

.0743 

.0797 

Table  16  shows  the  magnitude  of  the  projections  of  the 
estimated  smallest  right  singular  vectors  obtained  using 
inverse  iteration  (INV.  ITER.}  method  and  the  ICE,  onto  the 
smallest  singular  vector  computed  via  EVD. 

A  hypothesis  test  was  again  performed  on  these  results. 
Table  17  presents  the  results  obtained  for  the  hypothesis 
test . 


Table  17 :  RESULTS  FOR  THE  HYPOTHESIS  TEST  CONDUCTED 
ON  THE  COMPARISON  BETWEEN  THE  SV  INVERSE  ITERATION 
AND  THE  INCREMENTAL  CONDITION  ESTIMATOR. 


Test  Case 

Value  of  z 

Conclusion 

1 

-12.29 

Reject  Hq  =>  ^ 

(INV.  ITER.)  <  fi2  (ICE) 

2 

-4.74 

Reject  Hq  =>  fix  <  fi2 

3 

-3.76 

Reject  Ho  =>  fix  <  fi2 

4 

-  .  81 

Cannot  Reject  Ho 

36 


Tables  12  through  17  show  that  the  results  using  the 
incremental  condition  estimator  are  not  as  good  as  those  for 
the  SV  inverse  iteration.  Therefore,  the  inverse  iteration 
used  to  find  the  least  dominant  singular  vectors  is  preferred. 
The  source  code  implemented  to  approximate  the  smallest 
singular  value  and  its  corresponding  singular  vectors  via  ICE 
is  shown  in  Appendix  I. 
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IV.  COMPARISON  BETWEEN  THE  PERFORMANCE  OF  THE  RRQR  ALGORITHM 
ITERATED  FROM  DIMENSION  n  UNTIL  n-m+1  AND  FROM  DIMENSION  n 
UNTIL  n-r+1 

Prasad  [Ref.  5]  states  that  the  RRQR-based  algorithm  may 
be  used  to  estimate  the  DOA  information  when  the  algorithm  is 
iterated  from  n  until  n-m+1  rather  than  until  n-r+1  as  in 
Chan's  work  [Ref.  4].  The  relative  efficiency  of  the  RRQR 
algorithm  using  both  approaches  is  compared  by  running  1000 
trials.  The  scenario  is  m=2  fixed  sources  at  30°  and  32°  with 
n=10  sensors.  Therefore,  the  theoretical  noise -free 
autocorrelation  matrix  is .  of  size  10  x  10  and  has  rank  two 
(m=2)  .  The  near  rank  deficiency  is  (r=n-m)  8.  Three  SNR 

cases  are  tested  (-10,  0  and  10  dB)  for  each  one  of  the 

situations,  resulting  in  six  simulations. 

In  the  first  situation,  n-r+1  equals  three.  In  the 
second,  n-m+1  equals  nine.  It  is  intuitively  obvious  that  the 
second  one  is  much  faster  than  the  first,  as  it  will  be 
iterated  only  twice,  from  ten  to  nine.  On  the  other  hand,  the 
smaller  number  of  iterations,  the  less  probable  that  the  upper 
triangular  matrix  R  will  capture  the  near  rank  deficiency  of 
the  noise- free  autocorrelation  matrix  (R,)  . 

For  each  run,  the  largest  principal  angle  between  the 
signal  subspace  computed  via  the  eigenvector  decomposition  and 
the  approximated  signal  subspace  computed  via  the  RRQR 
algorithm  is  used  for  comparison.  The  same  is  done  for  the 
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noise  subspace.  A  vector  of  ones  is  expected  if  the  subspace 
generated  by  the  RRQR  algorithm  is  parallel  to  its  correspond¬ 
ing  subspace  generated  by  the  true  eigenvector.  In  the  case 
that  the  two  subspaces  are  perpendicular,  a  vector  of  zeros  is 
expected. 

Recall  that  the  cosines  of  the  principal  angles  between 
two  subspaces  F  and  G  are  defined  as  the  singular  values  of 
the  product  QFTQG  [Ref.  7]  ,  where  the  matrices  Qp  and  Qa  are  the 
orthonormal  matrices  obtained  from  F  and  G.  To  find  the  SVD 
decomposition  of  this  product,  the  inverse  cosine  of  the 
singular  values  is  taken  using  the  largest  angle.  The  largest 
angle  represents  a  measure  of  the  distance  between  the  two 
subspaces.  This  measure  is  used  for  noise  and  signal  subspa¬ 
ces  . 

Tables  18-20  present  means  and  standard  deviations 
obtained  for  the  largest  principal  angle  for  signal  and  noise 
subspaces  using  the  RRQR  algorithm  where  the  first  QR  decompo¬ 
sition,  in  step  0,  is  performed  with  pivoting.  This  computa¬ 
tion  is  done  for  SNR  -10,  0  and  10  dB. 
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Table  18:  COMPARISON  BETWEEN  THE  RRQR  ALGORITHM  ITERATED  FROM 
n  THROUGH  n-r+1  AND  FROM  n  THROUGH  n-m+1  FOR  SNR* -10  dB.  QR 
DECOMPOSITION  IN  STEP  0  PERFORMED  W/  PIVOTING,  1000  TRIALS. 


Iteration  from  n=l0 
until  n-r+l=3 
(#  of 

flops-1,301,732)  - 
Angle  in  degrees. 

Iteration  from  n=10 
until  n-m+l=9  (#  of 
flops=837, 515)  - 
Angle  in  degrees. 

Mean  (pc,) 

Std  Dev 

Mean  (pc2) 

Std  Dev 

Signal 

Subspace 

46.80 

13.65 

47.05 

14.11 

Noise 

Subspace 

46.80 

13.65 

47.05 

14.11 

Table  19  :  COMPARISON  BETWEEN  THE  RRQR  ALGORITHM  ITERATED  FROM 
n  THROUGH  n-r+1  AND  FROM  n  THROUGH  n-m+1  FOR  SNR-0  dB.  QR 
DECOMPOSITION  IN  STEP  0  PERFORMED  W/  PIVOTING,  1000  TRIALS. 


Iteration  from  n=10 
until  n-r+l=3 
(#.  of 

flops-1,301,732)  - 
Angle  in  degrees. 

Iteration  from  n=10 
until  n-m+l=9  (#  of 
flops=837, 515)  - 
Angle  in  degrees. 

Mean  ( pc, ) 

Std  Dev 

Mean  (pt2) 

Std  Dev 

Signal 

Subspace 

17.53 

5.71 

26.95 

12.55 

Noise 

Subspace 

17.53 

5.71 

26.95 

12.55 

Table  20:  COMPARISON  BETWEEN  THE  RRQR  ALGORITHM  ITERATED  FROM 
n  THROUGH  n-r+1  AND  FROM  n  THROUGH  n-m+1  FOR  SNR-10  dB.  QR 
DECOMPOSITION  IN  STEP  0  PERFORMED  W/  PIVOTING,  1000  TRIALS. 


Iteration  from  n=10 
until  n-r+l=3 
(#  of 

flops=l, 301, 732)  - 
Angle  in  degrees. 

Iteration  from  n=10 
until  n-m+l=9  {#  of 
f lops=837, 515)  - 
Angle  in  degrees. 

Mean  (pc,) 

Std  Dev 

Mean  (pt2) 

Std  Dev  | 

Signal 

Subspace 

2.24 

0.69 

4.85 

mm 

Noise 

Subspace 

2.24 

0.69 

4.85 

!■ 
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Tables  18-20  show  that  the  noise  and  signal  subspaces 
yield  identical  means  and  standard  deviations.  This  was  to  be 
expected  as  the  noise  and  signal  subspaces  contain  the  same 
information. 

Table  21  shows  the  results  obtained  for  the  hypothesis 
test  performed  to  compare  the  two  methods  of  iteration. 


Table  21:  RESULTS  FOR  THE  HYPOTHESIS  TEST  CONDUCTED  ON  THE 
COMPARISON  BETWEEN  THE  RRQR  ALG.  ITERATED  FROM  n  THROUGH  n-r-t-1 
AND  FROM  n  THROUGH  n-m+1.  QR  DEC.  IN  STEP  0  W/  PIVOTING 


Conclusion 

-10 

-  .40 

Cannot  reject  Hq 

0 

-21.6 

Reject  H o  =>  /ij  <  fi2 

10 

-19.9 

Reject  Hq  =>  Hi  <  n2 

The  results  show  that  for  SNR  0  and  10  dB,  the  null 
hypothesis  is  rejected.  Therefore,  the  larger  number  of 
iterations  yield  better  results.  However,  for  the  case  of  -10 
dB  the  results  are  meaningless  due  to  the  small  signal  to 
noise  ratio.  They  do  not  lead  to  detection  of  the  signal 
embedded  in  the  noisy  environment. 

Tables  22-24  present  means  and  standard  deviations 
obtained  for  the  largest  principal  angle  for  signal  and  noise 
subspaces  using  the  RRQR  algorithm  where  the  first  QR  decompo¬ 
sition,  in  step  0,  is  performed  without  pivoting.  An  hypothe¬ 
sis  test  is  not  necessary  in  this  case  due  to  the  clear 
difference  in  results  found  for  the  different  method  of 
iterations . 
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Table  22 :  COMPARISON  BETWEEN  THE  RRQR  ALGORITHM  ITERATED  FROM 
n  THROUGH  n-r+1  AND  FROM  n  THROUGH  n-m+1  FOR  SNR=-10  dB.  QR 
DECOMPOSITION  IN  STEP  0  PERFORMED  W/O  PIVOTING,  1000  TRIALS. 


Iteration  from  n=10 
until  n-r+l=3 
(#  of 

flops-1,301,732)  - 
Angle  in  degrees. 

Iteration  from  n=10 
until  n-m+l=9  (#  of 
flops=837, 515)  - 
Angle  in  degrees. 

Mean  ( ju, ) 

Std  Dev 

Mean  (^2) 

Signal 

Subspace 

46.76 

13.58 

59.59 

13.79 

Noise 

Subspace 

46.76 

13.58 

59.59 

13.79 

Table  23  :  COMPARISON  BETWEEN  THE  RRQR  ALGORITHM  ITERATED  FROM 
n  THROUGH  n-r+1  AND  FROM  n  THROUGH  n-m+1  FOR  SNR-0  dB.  QR 
DECOMPOSITION  IN  STEP  0  PERFORMED  W/O  PIVOTING,  1000  TRIALS. 


Iteration  from  n=10 
until  n-r+l=3 
(#  of 

flops-1, 301,732)  - 
Angle  in  degrees. 

Iteration  from  n=10 
until  n-m+l=9  (#  of 
flops=837, 515)  - 

Angle  in  degrees. 

Mean  ( ji, ) 

Std  Dev 

Mean  (fi2) 

Std  Dev 

Signal 

Subspace 

17.45 

5.70 

50.81 

23.62 

Noise 

Subspace 

17.45 

5.70 

50.81 

23.62 

Table  24:  COMPARISON  BETWEEN  THE  RRQR  ALGORITHM  ITERATED  FROM 
n  THROUGH  n-r+1  AND  FROM  n  THROUGH  n-m+1  FOR  SNR-10  dB.  QR 
DECOMPOSITION  IN  STEP  0  PERFORMED  W/O  PIVOTING,  1000  TRIALS. 


Iteration  from  n=10 
until  n-r+l=3 
{#  of 

flops-1,301,732)  - 
Angle  in  degrees. 

Iteration  from  n=10 
until  n-m+l=9  (#  of 
flops=837, 515)  - 
Angle  in  degrees. 

Mean  (/*,) 

Std  Dev 

Mean  (/x2) 

Std  Dev 

Signal 

Subspace 

2.21 

0.63 

14.47 

8.72 

Noise 

Subspace 

2.21 

0.63 

14.47 

8.72 
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Tables  22  to  24  show  that  better  performance  is  obtained 
when  using  a  larger  number  of  iterations.  Note  that  no 
difference  in  performance  is  found  using  a  QR  decomposition 
with  or  without  pivoting  in  step  0  of  the  RRQR  algorithm  in 
the  case  using  more  iterations.  This  result  is  not  true  for 
the  second  case.  This  is  clear  when  Tables  18-20  are  compared 
to  their  correspondent  Tables  22-24.  A  hypothesis  test  is  not 
needed  to  verify  this,  due  to  the  proximity  of  the  results. 
Therefore,  the  QR  decomposition  without  pivoting  is  preferred. 
This  method  is  less  computationally  intensive  when  the  origi¬ 
nal  algorithm  is  used,  performed  with  the  total  number  of 
iterations.  Note  that  the  option  using  only  two  iterations 
does  not  present  the  same  performance  regarding  the  pivoting. 
Therefore,  if  one  chooses  this  option,  care  should  be 
exercised  when  evaluating  the  pivoting  needs. 
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V.  COMPARING  THE  PERFORMANCE  OF  THE  MINIMUM  NORM  AND  MUSIC 
SPECTRAL  ESTIMATORS 

This  Chapter  investigated  the  computation  of  estimates  of 
the  direction  of  arrival  of  signals  obtained  using  the  RRQR. 
Two  high- resolution  techniques,  the  MUSIC  (Multiple  Signal 
Classification)  [Ref.  2]  and  the  Minimum  Norm  [Ref.  11]  are 
evaluated  to  verify  their  adequacy  when  used  with  the  RRQR 
algorithm  for  DOA  estimation. 

Rao  [Ref.  12]  shows  that  the  Mean  Square  Error  (MSE)  of 
the  Minimum  Norm  estimator  is  smaller  than  the  MSE  of  the 
MUSIC  estimator.  The  MUSIC  spectral  estimator  is  based  on 
the  orthogonality  principle.  Therefore,  the  principal 
eigenvectors  {v1#  v2,  .  .  . ,  vm}  span  the  same  subspace  as  the 
signal  vectors  {elf  e2,  .  .  . , em}  [Ref .  2].  Thus,  the  signal 
vectors  are  orthogonal  to  all  vectors  in  the  noise  subspace. 
The  power  density  corresponding  to  the  sources  DOA  information 
is  given  by: 


^ music n 

2  |eH(6>-Vj|2 


(43) 


where  Vj  is  the  singular  vectors  of  the  autocorrelation  matrix, 
m  is  the  number  of  signals,  n  is  the  number  of  sensors  or  the 
size  of  R,,  and  e(0)  is  the  mode  vector  as  defined  in  (4)  . 
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The  angle  6  is  varied  in  fine  steps  from  0  to  2ir.  When  it 
corresponds  to  one  of  the  source  DOA  angles,  e(0)  will  be 
equal  to  e^  i=l,2,..,m.  Since  the  sum  in  the  denominator  of 
(43)  is  over  m+1  through  n,  we  have  the  singular  vectors 
corresponding  to  the  noise  subspace.  By  the  orthogonality 
principle,  the  product  in  the  denominator  of  Pmusic  will  tend 
to  zero  and  Pmusic  will  tend  to  infinity  [Ref.  2]  .  The  result 
is  a  peak  at  the  source  DOA  angles. 

The  Minimum  Norm  estimator  is  based  on  the  estimation  of 
0,,  the  source  DOA  angles,  from  the  eigenstructure  of  the 
autocorrelation  matrix.  Suppose  we  have  a  vector  d  so  that 
d=  [dl,d2,  .  .  .  ,dj  ,  where  n  is  the  number  of  sensors  in  the 
array.  If  this  vector  has  the  property  that  xiHd=0, 
i=l,2,...,m,  where  m  is  the  number  of  sources  present  and  x 
the  data  snapshot  vector  at  instant  i,  then  we  can  find  a 
polynomial  D(z)  as 

n 

D(z)  diZ'i ,  (44) 
i=0 

so  that  the  zeros  of  the  polynomial  lie  at  the  elements  of  the 
mode  vector  corresponding  to  the  source  DOA  angles. 

The  polynomial  roots  corresponding  to  the  DOA  angle 
locations  lie  on  the  unit  circle  for  the  autocorrelation 
matrix  when  additive  noise  is  not  present.  However,  the  m 
roots  of  the  polynomial  D(z),  corresponding  to  the  m  sources. 
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lie  near  but  not  on  the  actual  unit  circle  when  the  estimated 


noise  free  autocorrelation  matrix  R,  is  used. 

The  Minimum  Norm  Spectral  Estimator  presents  the  following 
advantages : 


1.  The  estimates  of  0it  the  source  DOA  angles,  are  more 
accurate  even  at  relative  low  SNR  values  when  compared  with 
other  procedures . 

2.  The  n-m  extraneous  zeros  of  D(z)  tend  to  be  uniformly 
distributed  within  the  unit  circle  and  have  less  tendency 
for  false  sources. [Ref.  11] 


By  less  tendency  to  false  sources  we  mean  that  the  zeros 
of  D (z)  corresponding  to  noise  are  much  smaller  in  magnitude 
than  the  ones  corresponding  to  sources .  The  method  imposes 
dlf  the  first  element  of  vector  d,  to  be  equal  to  1  and 
requires  that  the  quantity  Q  in  Equation  (45)  be  minimum. 


k=l 


<4 12 


(45) 


The  effect  of  this  minimization  forces  the  extraneous  n-m 
zeros  to  be  uniformly  distributed  inside  the  unit  circle  [Ref. 
11]  . 

Next,  let  E,  be  the  matrix  constructed  with  the  signal 
subspace  generated  by  EVD  or  SVD.  We  partition  E,  as  follows: 


46 


(46) 


where  g=[eu,  e21,  .  .  . ,  eml] T  has  the  first  elements  of  the  signal 
subspace  eigenvectors  and  E,  is  the  matrix  E,  with  the  first 
row  deleted.  It  can  be  shown  that  in  order  to  satisfy  the 
desired  minimization  and  forcing  dt=l  [Ref.  11]  ,  we  have: 

1 

Esg*/(i-gHg) 

Once  the  vector  d,  representing  the  coefficients  of  the 
polynomial  D(z),  is  determined  via  (47),  the  roots  are 
computed.  There  are  m  roots  corresponding  to  the  sources  that 
present  a  large  magnitude  compared  to  the  ones  corresponding 
to  noise.  These  roots  reveal  the  desired  DOA  angles  in  their 
phase  angles. 

Results  for  both,  MUSIC  and  Minimum  Norm  spectral  estima¬ 
tors  are  shown  in  Figures  1-4,  for  two  fixed  sources  at  30° 
and  32°.  It  can  be  seen  from  these  examples  that  the  Minimum 
Norm  spectral  estimator  starts  to  resolve  the  two  sources  for 
a  SNR  of  5  dB,  while  MUSIC  starts  to  resolve  only  for  a  SNR  of 
7  dB.  This  agrees  with  the  results  shown  in  [Ref.  12] . 
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COMPARISON  OP*  MINORM  AMO  MUSIC  SPEC.  EST. .  SN«-  1 


ANGLE  IN  DEGREES 


s 


COMPARISON  OF  MINORM  AMO  MUSIC  SPEC.  EST. .  SNR  — 2 


ANGLE  IN  DEOREES 


Figure  1:  Comparison  between  the  DOA  estimated  by  MUSIC  and 
MINNORM  for  two  standing  sources  located  at  30°  and  32°  for 
SNR-1  and  2  dB. 
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COMPARISON  OR"  MINORM  AN  O  MUSIC  SPEC.  EST. .  SNR  — 3 


ANCLE  IN  DEGREES 


s 


Figure  2  s  Comparison  between  the  DOA  estimated  by  MUSIC  and 
MINNORM  for  two  standing  sources  located  at  30°  and  32°  for 
SNR* 3  and  4  dB . 
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ANGLE  INI  DEGREES 


Figure  3:  Comparison  between  the  DOA  estimated  by  MUSIC  and 
MINNORM  for  two  standing  sources  located  at  30°  and  32°  for 
SNR=5  and  6  dB. 


COMPARISON  OF  MlfSlORM  AMO  MUSIC  SPEC.  CST..  SNR  — 7 


Figure  4:  Comparison  between  the  DOA  estimated  by  MUSIC  and 
MINNORM  for  two  standing  sources  located  at  30°  and  32°  for 
SNR-7  dB. 


VI.  USING  THE  ADAPTIVE  RRQR  ALGORITHM  TO  TRACK  A  MOVING 
SIGNAL 

In  a  real  case,  we  are  interested  in  being  able  to  detect 
and  track  a  moving  source.  In  this  case,  the  moving  source 
information  was  sampled  every  t,  interval  to  provide  an  update 
of  the  source  information.  The  array  of  sensors  will  receive 
the  signals  from  the  sources  leading  to  the  data  vector 
x=  [x,,X2,  .  .  .  ,xj  .  From  vector  x,  we  compute  the 
autocorrelation  matrix,  R„,  from  (7)  and  form  the  noise-free 
autocorrelation  matrix,  Rs,  by  subtracting  the  noise  informa¬ 
tion,  a2 I,  Next  we  apply  the  RRQR  algorithm  for  each  update 
to  identify  the  signal  or  noise  subspaces.  Finally,  the 
Minimum  Norm  estimator  may  be  applied  which  leads  to  the 
identification  of  the  m  source  DOA  angles  from  the  n-m 
extraneous  noise  zeros.  This  algorithm  is  presented  in  this 
chapter. 

A.  THE  ALGORITHM 

In  order  to  track  the  DOA  of  a  moving  signal,  a  noise -free 
autocorrelation  matrix  must  first  be  generated.  Next,  we 
compute  a  first  RRQR  factorization  and  find  the  matrices  Q,  R 
and  II.  The  noise- free  autocorrelation  matrix,  R^R^-a2!,  may 
then  be  updated  for  each  one  of  the  next  signal  positions  in 
time,  according  to  Equation  (48)  below 
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(48) 


P  nev_  p  old 
Kxx  ~KXX 


+x-x  H 


new 


-X'XH 


old 


where  x  is  a  (n  x  1)  dimensional  vector  containing  the  input 
signals  at  each  one  of  the  n  array  sensors.  The  update  of  the 
noise- free  autocorrelation  matrix  is  achieved  using  a  moving 
window  where  the  number  of  snapshots  used  to  compute  the 
noise- free  autocorrelation  matrix  is  constant.  The  new  vector 
x  is  incorporated  in  the  autocorrelation  matrix  information 
for  each  snapshot  update,  while  the  oldest  snapshot 
information  is  discarded. 

A  possible  approach  in  updating  the  information  is  to  find 
the  updated  noise- free  autocorrelation  matrix  using  (48)  . 
Next  we  apply  a  QR  decomposition  followed  by  Chan's  [Ref.  4] 
pivoting  scheme,  i.e.,  a  complete  RRQR  algorithm. 

However,  it  is  possible  to  update  directly  the  existing  QR 
factorization  [Ref.  3]  for  each  new  time  sample.  Suppose  we 
want  to  add  a  rank- one  matrix  C=x-xH  to  the  matrix  R„old,  whose 
QR  factorization  is  known  as  QRlF.  The  new  matrix 
Rwnew=R«old+x •  xH  will  have  a  QR  factorization  Q1RiIIr,  where  x-x-1 
is  the  desired  rank- one  modification. 

The  rank- one  QR  factorization  update  may  be  computed  as 
follows : 

Rxxnew=  R]0told+x  •  xH  =  Q(R+QHx-xHn)lT=  QiRjIT 

Let  w=Qh-x.  Complex  Givens  rotations  can  be  used  to  zero 
out  all  elements  of  w  except  for  the  first  component,  leading 
to 


52 


(49) 


J1H.  .  .Ja-iW  =  [a  0  0  ...  0]  T 
If  the  same  Givens  rotations  are  applied  to  R,  it  can  be  shown 
that  H=JjH.  .  .  Jn_tHR  is  upper  Hessenberg.  Consequently,  we  have 
( J,H.  .  . Jn.!H)  (R„ol<1+wxH)  =H+ [a  0  ...  0]t-xh=H,,  also  upper 

Hessenberg. [Ref .  7] 

Next,  we  use  complex  Givens  rotations  to  compute 
G,”.  .  .Gn.1HH1=R1,  where  R,  is  upper  triangular.  Combining  all  of 
the  above,  we  have  the  desired  new  QR  factorization 
Rxxnew=Rxxold+x  •  xh=Q,  •  Rr IF ,  where 

Q1=Q‘Jn.1 .  .  .  Jx ‘G1 .  .  .  Gn_x .  (50) 

The  reader  should  refer  to  Golub  [Ref.  7]  for  additional 
detail . 

We  apply  two  successive  rank- one  modifications  to  update 
the  noise- free  autocorrelation  matrix  for  the  current  time 
sample,  finding  a  new  set  of  matrices  Q,  R  and  II.  The  first 
rank- one  modification  incorporates  the  new  data  vector  to  the 
autocorrelation  matrix.  The  second  accounts  for  removing  the 
old  information.  Then,  we  apply  Chan's  [Ref.  4]  pivoting 
scheme  to  insure  a  correct  estimation  of  the  noise  and/or 
signal  subspaces. 

Next,  we  identify  the  signal  and  noise  subspaces.  The 
first  m  columns  of  Q  constitute  the  signal  subspace,  where  m 
is  the  number  of  sources  present.  The  last  n-m=r  columns 
constitute  the  noise  subspace,  where  n  is  the  number  of 
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sensors.  Note  that  it  is  usually  more  efficient  to  use  the 
signal  subspace,  rather  than  the  noise  subspace  as  it  is 
smaller.  Finally,  the  Minimum  Norm  algorithm  is  applied  to 
estimate  the  source  locations. 

A  total  of  m  out  of  the  n  roots  of  the  polynomial  whose 
coefficients  form  the  vector  "d"  in  the  MINORM  algorithm 
corresponds  to  the  source  locations.  It  is  expected  that  the 
m  roots  corresponding  to  the  sources  lie  near  the  unit  circle 
and  the  remaining  ones  have  smaller  magnitudes.  Some  sort  of 
filtering  must  be  applied  to  separate  the  m  source  zeros  and 
the  remaining  n-m  extraneous  zeros.  Filtering  may  be  achieved 
either  by  sorting  the  expected  range  of  source  angles  and/or 
by  sorting  the  magnitude.  Thus,  the  algorithm  corresponding 
to  the  adaptive  case  becomes  fRef .  3] : 


1)  Initialization 

Rs=X-Xh-w<t2I  (note  that  the  noise- free  autocorrelation 
matrix  is  unnormalized)  where  X  is  defined  in  Equation 
(7) ,  Chapter  1  and  w  is  the  number  of  snapshots  used  to 
compute  the  correlation  matrix. 


RRQR  factorization  of  R,  (RJI=QR) 

2)  Start  updating.  For  k*l  to  number  of  updates  do: 

a)  R,(k+1)  =R,(k)  +x(k+l)  -xH(k+l)  -x(k+l-w)  -xH(k+l-w)  . 

b)  Update  the  above  QR  factorization  applying  two  succes¬ 
sive  rank-one  modifications  to  R,(k),  leading  to: 

R, (k+1)  =0^^  (Alternatively  we  may  find  a  new  complete 
RRQR  decomposition  of  the  updated  noise- free 
autocorrelation  matrix.  In  this  case,  the  next  two  steps 
should  be  skipped) . 
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c)  Apply  Chan's  Pivoting  scheme  (steps  l  through  9  of  the 

RRQR  algorithm)  to  II,  Q[  and  R, ,  finding  IIj,  Q2  and  R2. 

d)  Let  11=1X2 ,  Q=Q2  and  R=R2. 

e)  Identify  the  signal  or  noise  subspace  (Option  to  use  a 

refinement  as  explained  in  Section  B  below) . 

f)  Apply  the  Minimum  Norm  to  find  the  estimated  source 

angles . 

g)  Filter  and  store  the  source  angles. 

The  source  code  implemented  to  generate  the 
autocorrelation  matrix  is  shown  in  Appendices  K  and  L.  The 
source  code  implemented  to  generate  the  adaptive  algorithm  is 
shown  in  Appendix  D  and  the  code  corresponding  to  the  rank -one 
modification  shown  in  Appendix  E.  Appendix  F  presents  the 
Minimum  Norm  algorithm.  Finally,  Appendix  H  presents  the 
Minimum  Norm  identification  procedure  needed  to  isolate  the 
signal  source  locations. 

The  adaptive  algorithm  presented  above  is  an  alternative 
to  the  identification  of  the  signal/noise  subspaces  via  SVD  or 
EVD  decomposition.  Note  that  steps  4  to  6  in  Chan's  algorithm 
are  only  needed  when  the  element  of  maximum  magnitude  of  the 
smallest  right  singular  vector  is  not  in  the  last  position. 
Thus,  additional  reduction  in  computation  time  is  obtained 
when  no  pivoting  is  needed. 

Next,  we  investigate  how  often  pivoting  is  needed  in 
Chan's  algorithm.  To  test  that  effect,  we  ran  two  test  cases 
with  two  sources  (m=2)  .  One  of  the  sources  is  fixed,  the 
other  is  time  varying.  Ten  sensors  are  used  to  compute  the 
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correlation  matrix  with  100  snapshots  used  to  form  the  noise- 
free  autocorrelation  matrix  and  200  updates  are  computed.  The 
value  of  the  SNR  varies  to  try  to  identify  correlation  between 
the  SNR  and  the  need  for  pivoting  during  the  RRQR  decomposi¬ 
tion  (step  4  of  Chan's  algorithm).  In  our  test  cases,  if  the 
Chan's  pivoting  scheme  were  applied  blindly,  steps  4  through 
9  would  be  executed  for  a  total  of  1600  times.  (Because  n- 
r+l=10-8+l=3  and  the  algorithm  is  iterated  from  n  to  n-r+1,  a 
total  of  eight  times  for  each  one  of  the  updates  is  needed. 
Since  we  have  200  updates  in  our  test  case,  we  get 
8*200=1600) . 

First,  we  update  Q  and  R  directly  using  two  successive 
rank-one  modifications.  Next,  we  update  the  noise-free 
autocorrelation  matrix  using  Equation  (48)  and  perform  a  new 
complete  RRQR  decomposition  (steps  0  through  9) .  Figures  5 
and  6  show  the  percentage  of  times  that  a  pivoting  is  needed 
out  of  the  1600  tests  as  a  function  of  the  SNR  (dB)  of  the 
source . 

Figure  5  shows  that  the  percentage  of  pivoting  steps 
needed  lies  between  18%  and  32%  when  two  successive  rank-one 
modifications  are  used  to  update  the  noise -free  correlation 
matrix.  This  means  that  no  pivoting  is  needed  for  every 
iteration. 

For  the  second  test  case,  where  a  complete  RRQR  algorithm 
is  applied  in  the  updated  noise- free  autocorrelation  matrix, 
the  percentage  of  pivoting  steps  needed  remains  between  70% 
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PERCENTAGE  OF  PIVOTING  NEEDED  FOR  200  UPDATES 


Figure  5:  Percentage  of  times  that  a  pivoting  is  needed  out  of 
1600  iterations  on  the  adaptive  algorithm,  using  two  rank-one 
modifications. 


PERCENTAGE  OP"  PIVOTI  MG  NEEDED  EQR  200  UPDATES 


Figure  6:  Percentage  of  times  that  pivoting  is  needed  out  of 
1600  iterations  on  the  adaptive  alg.,  updating  the  noise-free 
correlation  matrix  and  applying  a  complete  RRQR  algorithm. 

and  88%.  The  two  successive  rank-one  modifications  are  more 
computationally  expensive  than  the  complete  RRQR  algorithm. 
On  the  other  hand,  they  require  less  pivoting.  In  Section  D 
below,  we  analyze  the  implications  of  this  on  the  processing 
time  for  the  two  approaches.  Furthermore,  Figures  5  and  6 
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indicate  that  there  is  no  correlation  between  the  magnitude  of 
the  SNR  and  the  need  for  pivoting. 


B.  REFINEMENT  ON  THE  RESULTS  OF  THE  ADAPTIVE  CASE 

This  section  presents  an  improvement  to  the  signal 
subspace  estimation  procedure  which  is  applied  to  the  adaptive 
algorithm.  The  basic  idea  behind  the  improvement  is  based  on 
the  fact  that  the  noise- free  autocorrelation  matrix,  R,=QRllT 
is  Hermitian,  and  therefore  RS=R,H.  So,  R,H=IIRHQH.  Recall  that 
the  signal  and  noise  subspaces,  Qs  and  Qn  are  contained  in  the 
matrix  Q.  Thus, 


RsO  ,=IIR* 


Qs 

Qn 


H 


H  P.°| 


=  hr11 


Is  0 
0  0 


n 


is"  < 

iff"  ff"i 
l«12  KZ2 


Is  O' 
0  0 


=  n 


l*ii  o 
1*12  0 


(51) 


and  therefore 


RSQS  =  (nx  n2j 


1*11  0 
1*12  0 


=  n^+n^. 


(52) 


The  matrix  resulting  by  the  product  R,Q,  may  be  viewed  as 
an  one- step  subspace  iteration  applied  to  the  signal  subspace 
Q,.  This  iteration  scheme  improves  the  results  obtained  as 
shown  in  the  next  section.  The  drawback  is  that  the  resulting 
iterated  signal  subspace  R,Q,  is  no  longer  orthonormal.  An 
additional  orthonormalization  step  needs  to  be  applied  to  the 
iterated  signal  subspace  in  order  to  use  the  Minimum  Norm 
algorithm.  Note  that  no  reorthonormalization  is  needed  when 
the  MUSIC  estimator  is  used. 
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C.  SIMULATION  RESULTS 


We  now  present  the  results  generated  by  the  simulation 
carried  out  for  SNRs  of  6,  10  and  20  dB.  The  noise  is  assumed 
to  be  zero -mean  Gaussian  and  uncorrelated  from  sensor  to 
sensor.  We  consider  the  case  of  two  sources  impinging  on  the 
ten- element  array.  The  first  source  is  assumed  to  be  fixed  at 
a  normalized  angle  #[=40°,  the  second  source  location  02  is 
linear  time-varying.  Movement  starts  at  30°  and  stops  at 
22.5°,  after  100  snapshots  are  used  to  form  the  correlation 
matrix  and  200  updates  are  used  to  simulate  the  movement. 

Figures  7,  8  and  9  present  the  estimated  DOA  information 
obtained  using  the  Eigenvector  decomposition  (EVD) ,  the 
initial  RRQR  approximation,  and  the  "refined"  RRQR  algorithm. 
The  EVD  decomposition  is  used  for  convenience  instead  of  the 
SVD  decomposition,  as  both  span  the  same  subspace.  The 
results  obtained  for  the  refined  RRQR  technique  are  nearly 
identical  to  those  obtained  using  the  EVD  technique.  Table  25 
below  shows  means  and  standard  deviations  for  the  magnitude  of 
the  difference  between  the  RRQR  approximation  for  the  signal 
DOA  angle  in  degrees  with/without  refinement  and  the  DOA 
generated  by  the  EVD  decomposition. 
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LOCATION  OF  MOVING  SOURCE  -  SNR=6  dB 
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Figure  7:  Position  of  the  time-varying  source  for  each  update, 
SNR  equal  to  6  dB. 


60 


NB  OF  UPDATES 


LOCATION  OF  MOVING  SOURCE  -  SNR-1 0  dB 


Nouvooi  aoanos 


Figure  8:  Position  of  the  time- varying  source  for  each  update, 
SNR  equal  to  10  dB. 
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NB  OF  UPDATES 


LOCATION  OF  MOVING  SOURCE  -  SNR-20  dB 


NoiivDoi  aoanos 


Figure  9:  Position  of  the  time-varying  source  for  each  update, 
SNR  equal  to  20  dB . 
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OF  UPDATES 


Table  25:  MEAN  AND  STD.  DEV.  FOR  THE  MAGNITUDE  OF  THE  DIFFER¬ 
ENCE  BETWEEN  THE  RRQR  APPROX.  FOR  THE  SIGNAL  DOA  ANGLE  IN  DEG. 
FOR  THE  MOVING  SOURCE  WITH/WITHOUT  REFINEMENT  AND  THE  EVD. 


SNR=6  dB 

SNR=10  dB 

SNR=20  dB 

Average  w/o 

refinement 

(degrees) 

1.2346 

.4530 

.0525 

Standard  Devi¬ 
ation  w/o  re- 
f inement  ( de - 
grees) 

.7569 

.2853 

.0302 

Average  w/ 
refinement 
(degrees) 

.0542 

.0079 

8.6930  E-5 

Standard  Devi¬ 
ation  w/  re¬ 
finement  (de¬ 
grees) 

.0323 

.0061  . 

6.2626  E-5 

The  results  are  excellent  even  for  the  approximation  without 
refinement.  The  refinement  improvement  becomes  better  as  the 
S^  increases. 

Results  found  for  the  largest  principal  angle  between  the 
projection  of  the  signal  subspace  found  via  RRQR  and  the  true 
signal  subspace  found  via  EVD  are  shown  next.  Figures  10,  11 
and  12  depict  the  results  found  for  SNR  values  equal  to  6,  10 
and  20  dB. 
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PRINC.  ANGLE  FOR  THE  SIG.  SBSP-SNR.6  dB 


Figure  10 :  Largest  principal  angle  between  the  signal  subspace 
found  via  RRQR  and  the  true  signal  subspace  found  via  EVD  for 
SNR* 6  dB. 


PRINC.  ANGLE  FOR  THE  SIG.  SBSP-SNR-10  dB 
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UPDATE  NUMBER  SOLID  *W/0  REFINEMENT 
UPDATE  NUMBER  OASHED-W/  REFINEMENT 


Figure  11:  Largest  principal  angle  between  the  signal  subspace 
found  via  RRQR  and  the  true  signal  subspace  found  via  EVD  for 
SNR-10  dB. 


64 
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PRINC.  ANGLE  FOR  THE  SIC.  SBSP-SNR-20  <JB 


UPDATE  NUMBER 


SOLID  -W/O  REFINEMENT 
DASHED-W/  REFINEMENT 


Figure  12 :  Largest  principal  angle  between  the  signal  subspace 
found  via  RRQR  and  the  true  signal  subspace  found  via  EVD  for 
SNR=20  dB. 

As  can  be  seen,  the  results  for  the  cases  with  refinement  are 
much  better  than  those  without  refinement.  Table  26  shows  the 
means  and  standard  deviations  obtained  for  the  largest 
principal  angle  between  the  RRQR  and  EVD  for  both  cases,  with 
and  without  refinement. 
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Table  26:  MEAN  AND  STANDARD  DEVIATION  FOR  THE  LARGEST 
PRINCIPAL  ANGLE  IN  DEGREES  BETWEEN  THE  RRQR  APPROXIMATION  AND 
THE  EVD  FOR  THE  SIGNAL  SUBSPACE  WITH/WITHOUT  REFINEMENT. 


1  _ 

SNR =6  dB 

SNR=10  dB 

SNR=20  dB 

Him 

.81707 

3.4087 

.3442 

Standard  Devi¬ 
ation  w/o  re¬ 
finement 
(degrees) 

1.5944 

.6502 

.0595 

1.2878 

.2252 

.0023 

Standard  Devi¬ 
ation  w/  re- 
f inement 
(degrees) 

.4368 

.0737 

6.9273-4 

Last,  Figures  13,  14  and  15  present  the  behavior  of  the 
RRQR  approximation  for  the  estimation  of  the  DOA  of  the  second 
source  that  remained  constant  at  40°.  Table  27  presents  mean 
and  standard  deviation  values  for  the  magnitude  of  the 
difference  between  the  DOA  estimated  by  the  RRQR  with/without 
refinement  and  the  EVD.  As  can  be  seen,  the  RRQR  results 
agree  closely  with  EVD  results. 
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SIGNAL  SUBSPACE  -  SNR-6  dB 


nouvdoi  aoanos 


Figure  13:  DOA  in  degrees  for  the  fixed  source  for 
update,  SNR  equal  to  6  dB. 


each 
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NB  OF  UPDATES 


SIGNAL  SUBSPACE  -  SNR-10  dB 


nouvooi  aoanos 


Figure  14:  DOA  in  degrees  for  the  fixed  source  for  each 
update,  SNR  equal  to  10  dB. 
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NB  OF  UPDATES 


SIGNAL  SUBSPACE  -  SNR-20  dB 


noiivooi  aoanos 


Figure  15:  OOA  in  degrees  of  the  fixed  source  for  each  update, 
SNR  equal  to  20  dB. 
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NB  OF  UPDATES 


Table  27:  MEAN  AND  STD.  DEV.  FOR  THE  MAGNITUDE  OF  THE  DIFFER¬ 
ENCE  BETWEEN  THE  RRQR  APPROX.  FOR  THE  SIG.  DOA  ANGLE  IN  DEG. 
WITH/WITHOUT  REFINEMENT  AND  THE  EVD  FOR  THE  FIXED  SOURCE. 


SNR=6  dB 

SNR=10  dB 

SNR=20  dB 

Average  w/o  re¬ 
finement  (degrees) 

.  8970 

.3623 

.  0418 

Standard  Deviation 
w/o  refinement 
(degrees) 

.4603 

.1619 

.  0184 

Average  w/  re¬ 
finement  (degrees) 

.0226 

.0049 

7 . 0319E- 5 

Standard  Deviation 
w/  refinement 
(degrees) 

.0201 

.0047 

6 . 4206E- 5 

D.  COMPUTATIONAL  EFFICIENCY  OF  THE  RRQR  APPROXIMATION 

This  section  presents  a  basic  estimation  of  the  number  of 
floating-point -operations  (flops)  needed  to  compute  the  RRQR 
approximation.  The  n- dimensional  noise -free  autocorrelation 
matrix  is  square  and  is  not  considered  complex  valued  for  flop 
computation.  This  fact  will  be  taken  into  account  later.  The 
number  of  flops  necessary  to  perform  one  SVD  decomposition  is 
6n3  [Ref.  4]  .  Every  time  a  new  sample  arrives,  we  need  an 
0(n3)  operation  to  recompute  the  SVD  [Ref.  13]. 

The  RRQR  algorithm  is  composed  of  three  distinct  parts. 
The  computation  of  the  initial  QR  factorization  without 
pivoting,  computed  only  once  in  the  beginning  of  the  algo¬ 
rithm;  the  computation  of  the  least  dominant  right  singular 
vector  v  of  Rn  by  inverse  iteration,  at  each  iteration;  and 
the  new  QR  factorization  of  RnP,  also  at  each  iteration. 


The  first  part  takes  2n3/3  flops  using  the  Householder 
algorithm  if  we  do  not  need  to  accumulate  Q.  If  Q  is  needed, 
as  in  our  case,  it  takes  4n3/3  flops  [Ref.  7]  .  The  Modified 
Gram-Schmidt  algorithm  is  preferred,  since  it  takes  only  n3 
flops  when  the  accumulation  of  Q  is  performed  [Ref.  7] . 

Ignoring  lower  order  terms,  the  second  part  of  the  RRQR 
factorization  takes  In2r  flops  to  be  iterated  [Ref.  4],  where 
r  is  the  noise- free  autocorrelation  matrix  rank -deficiency  and 
I  is  the  number  of  iterations  used  in  the  inverse  iteration 
method.  Our  case  uses  1=3,  therefore,  the  second  part  takes 
3n2r  flops  to  be  performed. 

In  the  third  part  of  the  RRQR  algorithm,  2n2r  flops  are 
needed  when  Givens  rotations  are  used  [Ref.  4]  .  However,  note 
that  not  all  elements  below  the  main  diagonal  of  the  matrix 
RUP  needs  to  b"  ihilated  because  they  are  already  zero. 
Therefore,  a  conditional  "IF”  statement  may  be  used  to  verify 
if  the  element  is  already  zero  to  save  additional  flops. 
Thus,  the  RRQR  algorithm  totals  n3+5n2r  flops  at  most. 

Following  the  first  RRQR  decomposition,  we  have  two 
options.  The  first  option  updates  the  autocorrelation  matrix, 
as  in  Equation  (48)  and  performs  a  new  QR  decomposition.  This 
update  is  followed  by  the  pivoting  scheme,  as  in  steps  0 
through  9  of  Chan's  algorithm.  The  second  option  updates  the 
already  obtained  QR  decomposition  directly,  using  two 
successive  rank-one  modifications.  A  new  Q,  and  R,  is  found 
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and  the  Chan's  pivoting  scheme  is  applied.  An  analysis  of  the 
number  of  flops  necessary  to  perform  the  two  options  is  made. 

Recall  that  in  the  adaptive  case,  a  double  update  must  be 
performed  for  each  sample.  One  update  adds  the  new  sample  and 
the  other  subtracts  the  old  one.  If  the  rank- one  modification 
option  is  used,  the  updates  take  13n2  flops  each  for  each 
snapshot,  totaling  26n2  flops  [Ref.  7].  This  does  not  include 
the  number  of  flops  necessary  to  proceed  the  pivoting  scheme. 
If  the  complete  RRQR  factorization  option  is  used,  the  QR 
decomposition  using  the  Modified  Gram- Schmidt  algorithm  takes 
n3  flops.  Note  that  in  such  a  case,  the  autocorrelation 
matrix  must  be  updated  as  Rnew=Rold+ (x-xH) new-  (x-xH)old,  which  takes 
4n2  flops  (n2  for  each  multiplication  and  n2  for  each 
addition)  .  The  total  is  n3+4n2  for  each  snapshot. 

Comparing  n3+4n2  with  26n2,  we  see  that  to  update  the 
noise- free  autocorrelation  matrix  (finding  R^  as  in  Equation 
(48))  and  to  take  its  QR  decomposition  is  more  economical  than 
to  perform  two  rank- one  modifications  for  0  <  n  <  22. 
Therefore,  a  complete  RRQR  algorithm  including  a  Modified 
Gram- Schmidt  QR  decomposition  and  the  pivoting  scheme  is 
preferable  when  compared  to  the  two  rank-one  modification 
updates,  for  a  number  of  sensors  smaller  than  22.  This  method 
does  not  take  into  consideration  the  difference  of  pivoting 
schemes  needed  after  the  update.  Table  28  presents  a 
comparison  for  the  number  of  flops  needed  for  both  processes. 
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Table  28:  COMPARISON  OF  THE  NUMBER  OF  FLOPS  NEEDED  TO  THE 
SNAPSHOT  UPDATE  OF  THE  ADAPTIVE  ALGORITHM  VIA  A  COMPLETE  RRQR 
FACTORIZATION  AND  TWO  RANK -ONE  MODIFICATIONS. 


Snapshot  update 
via  complete 
RRQR  algorithm 

Number  of 
flops  need¬ 
ed 

Snapshot  update 
via  two  rank- 
one  modif. 

Number  of 
flops  need¬ 
ed 

Computation  of 

^new 

4n2 

13n" 

QR  Factoriz.  of 

Rncw 

nJ 

Second  rank- one 
modif . 

13ni_ 

Smallest  right 
SV 

3rPr 

Smallest  right 

SV 

3n'r 

iiii  *1? 

—  2a*r 

2x7t 

Total 

Total 

As  seen  earlier,  using  two  successive  rank- one 
modifications  leads  to  pivoting  for  at  most  32%  of  the  1600 
iterations  used  to  perform  the  200  updates.  This  is  compared 
to  the  maximum  of  88%  pivoting  when  updating  the  noise- free 
autocorrelation  matrix  and  applying  a  complete  RRQR  algorithm. 
These  numbers  are  used  to  compare  the  two  approaches. 

For  the  third  part  of  the  RRQR  algorithm,  we  need  2n2r 
flops.  Using  the  first  approach,  two  rank- one  modifications 
over  Q  and  R  requires  a  total  of  2  6n2+ (3+0 . 32x2)  n2r  flops  for 
each  update.  In  the  second  approach,  the  updating  of  the 
noise- free  autocorrelation  matrix  and  the  application  of  a 
complete  RRQR  algorithm  requires  n3+4n2+ (3+0 . 88x2)  n2r  flops  for 
each  update.  Figure  16  depicts  the  results  showing  the 
regions  when  the  complete  RRQR  or  the  update  approaches  might 
be  preferred. 
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Regions  of  best  use  for  the  RRQR  and  Direct  update  approaches 
- , - ! - 1 - r 


0  5  tO  15  20 


RANK  DEFICIENCY  (r) 

Figure  16:  Graph  showing  the  regions  where  the  RRQR  or  the  two 
rank-one  modifications  approach  is  preferred* 

The  MINNORM  algorithm  takes  roughly  a  number  of  flops 
equals  to  2nr  plus  the  flops  necessary  to  compute  the 
polynomial  roots.  The  root  finding  routine  is  iterative. 
Thus,  it  is  impossible  to  obtain  a  specific  expression  to 
evaluate  the  number  of  flops  necessary  to  run  it.  However, 
for  the  scenario  simulated  in  the  previous  section,  we  were 
able  to  compute  the  mean  and  standard  deviation  of  the  number 
of  flops  spent  by  the  MATLAB™  software  to  find  the  ten  roots 
for  each  one  of  the  200  updates. 

To  evaluate  the  sensitivity  of  the  root  finding  routine  to 
the  SNR,  we  tested  for  SNR=-100,  6  and  100  dB.  The  results 
are  shown  in  Table  29. 
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Table  29:  MEAN  AND  STD.  DEV.  FOR  THE  NUMBER  OF  FLOPS  NECESSARY 
TO  FIND  TEN  ROOTS  OF  THE  POLYNOMIAL  FORMED  BY  THE  MINIMUM  NORM 
ALGORITHM.  SNR- -100,  6,  100  DB. 


SNR (dB) 

Mean  of 
#  flops 

Std. 
Dev.  # 
flops 

-100 

63159 

3042 

6 

61007 

2168 

100 

64420 

148 

According  to  Table  29,  the  number  of  flops  necessary  to 
the  root  finding  routine  presents  a  smaller  standard  deviation 
for  SNR  equals  to  100  dB.  This  happens  because  at  such  a  high 
SNR,  the  zeros  are  at  more  or  less  fixed  locations.  Thus,  one 
can  expect  about  the  same  number  of  iterations  needed  to 
identify  them.  The  study  of  the  polynomial  root  algorithm 
sensitiveness  to  the  number  of  sources  at  different  locations 
deserves  further  research  and  is  out  of  the  scope  of  this 
work. 

The  REFINEMENT  algorithm  spends  a  total  of  3n3+ (l-2r)  n2-rn 
flops,  including  the  orthonormalization  necessary  to  be  used 
in  conjunction  with  the  MINNORM.  When  the  problem  is  located 
under  the  line  depicted  in  Figure  16,  the  updating  via  a 
complete  RRQR  algorithm  is  preferred  in  the  most  probable 
case.  Totalizing,  the  RRQR,  the  MINNORM  and  the  REFINEMENT 
algorithms  take  a  total  number  of  flops  of  4n3+ (2 . 76r+5)  n2+rn 
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for  each  update.  This  value  was  developed  for  a  real  valued 
noise- free  autocorrelation  matrix.  To  cope  with  future  growth 
of  the  program,  we  might  multiply  the  number  of  flops  by  a  4 
to  1  factor  when  implementing  it  operationally  and  by  a 
roughly  4  to  1  factor  to  estimate  the  case  of  a  complex 
matrix. 

E.  INTEGRATING  THE  RRQR  ALGORITHM  INTO  A  REAL-TIME  CASE 

Suppose  we  have  a  signal  being  received  by  a  linear  phased 
array  on  the  earth  surface  from  a  moving  object  located  at  10 
NM  from  the  sensors  at  an  angle  of  80°  from  the  sensors 
vertical  (see  Figure  17) . 


80° 


n  sensors 


m  <  n 


m  sources 


Passive  Linear 

Equispaced 

Array. 


Figure  17 :  Situation  of  a  signal  being  received  by  linear 
phased  array  sensors. 
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Also  suppose  that  this  signal  is  moving  at  an  horizontal 
velocity  v.  Based  on  the  simulations  run  in  this  thesis,  the 
moving  signal  changed  5°  in  its  DOA  with  respect  to  the  array 
vertical.  This  was  done  in  200  snapshots.  Each  snapshot  is 
taken  at  an  interval  of  ts  seconds.  Therefore,  the  angular 
speed  of  this  moving  signal  is  defined  as 


180  id 
200 . ts  S 


The  signal  velocity  is 


(53) 


v= 


w-R 

COS (0) 


(54) 


where  R  is  the  distance  between  the  signal  and  the  array. 
Using  (53)  and  (54)  leads  to: 


£  _ _ 5  "it  'R _ 

s  36000‘VCOS (0) 


(55) 


Assuming  a  signal  is  moving  at  sound  speed  as  in  Figure  18, 
0=80 ° ,  R=10  NM,  a  t,  of  136.87  ms  would  be  needed. 

For  the  scenario  simulated  in  the  previous  section  where 
n=10  and  r=8,  the  number  of  flops  needed  for  each  update  would 
be  16x6788=108.6  Kflops  per  update.  This  value  is  equivalent 
to  108,600/136.87  ms  *»  790  Kflops  per  second  of  microprocessor 
computing  power.  Assuming  a  reasonable  computing  power  of  one 
flop  per  clock  pulse  at  32  bits,  a  microprocessor  would  have 
to  operate  at  approximately  790  KHz.  Current  commercial 
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microprocessors  sold  with  a  clock  of  50  MHz  would  be  able  to 
implement  this  algorithm  in  real-time.  Note  that  there  is 
enough  room  to  accommodate  the  root  finding  procedure 
neglected  in  the  MINNORM  algorithm  flops  computation. 


78 


VII. 


CONCLUSIONS 


We  have  presented  a  fast  algorithm  to  isolate  signal  and 
noise  subspaces  without  performing  the  eigenvector 
decomposition  of  the  autocorrelation  matrix.  This  method  is 
called  the  Rank -Revealing  QR  factorization. 

We  have  investigated  the  possibility  of  using  the  least 
dominant  eigenvector  instead  of  using  the  least  dominant 
singular  vector  in  step  two  of  the  RRQR  algorithm. 
Simulations  have  shown  that  the  minimum  singular  vector 
generated  by  the  inverse  iteration  method  gives  better  results 
than  those  obtained  with  the  approximate  smallest  eigenvector 
generated  by  the  same  method.  The  Incremental  Condition 
Estimator  algorithm  for  finding  an  approximation  for  the 
smallest  singular  vector  has  also  been  tested.  Simulations 
have  demonstrated  that  the  inverse  iteration  again  yields 
better  results. 

The  possibility  of  using  a  faster  RRQR  algorithm  than  the 
original  algorithm  with  fewer  iterations  has  been 
investigated.  Simulations  have  shown  that  reducing  the  number 
of  iterations  worsens  the  signal/noise  subspace  estimations. 
Therefore,  the  original  algorithm  is  preferred. 

Two  spectral  estimators  have  been  tested  to  be  used  with 
the  RRQR  algorithm,  the  MUSIC  and  the  Minimum  Norm.  A  limited 
number  of  simulations  have  indicated  that  the  Minimum  Norm 
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resolved  two  signals  for  a  lower  SNR  than  the  MUSIC  method,  in 
agreement  with  Rao  [Ref.  12]  . 

A  relatively  inexpensive  computational  refinement 
algorithm  has  been  presented  for  the  estimation  of  the  RRQR 
signal  subspace.  Simulations  have  shown  that  improvements  at 
least  as  high  as  a  factor  of  20  are  possible  to  be  obtained 
for  the  signal  angle  of  arrival,  as  compared  with  the  original 
RRQR-based  DOA  results.  Note  that  this  refinement  improvement 
becomes  better  as  the  SNR  increases . 

An  adaptive  RRQR-based  algorithm  has  been  introduced  to 
track  the  DOA  of  moving  signals.  Two  options  have  been 
evaluated  to  compute  correlation  updates.  The  more  adequate 
option  may  be  determined  depending  upon  the  particular  problem 
set  up,  e.g.,  the  noise-free  autocorrelation  matrix  rank 
deficiency  (r)  and  the  number  of  sensors  (n) . 

An  evaluation  of  the  number  of  flops  required  by  the 
adaptive  algorithm  to  find  the  DOA  for  two  sources  present  has 
been  determined.  The  results  have  shown  the  feasibility  of 
the  algorithm  to  solve  a  real-time  problem. 
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Appendix  A 


% 

%  This  function  implements  Chan's  RRQR  algorithm. 

%  Milton  P.  Ferreira,  Sep/1992. 

% 

%  Input  parameters: 

% 

%  r=  matrix  to  which  we  want  to  apply  the  RRQR  factorization 
%  (noise- free  autocorrelation  matrix. 

%  rl=  rank  deficiency  of  the  matrix  "r". 

% 

%  Output  parameters: 

% 

%  Q=  orthonormal  matrix  resulting  from  the  RRQR  factorization 
%  of  "r" . 

%  R-  upper  triangular  matrix  resulting  from  the  RRQR 
%  factorization  of  "r" . 

%  e=  permutation  matrix  (PI)  resulting  from  the  RRQR 
%  factorization  of  "r" . 

%  nsbsp=  noise  subspace. 

%  ssbsp=  signal  subspace. 

% 

function  [Q, R, e, nsbsp, ssbsp] =rrqr (r, rl) 
n=size (r) ;n=n(l, 1) ; 
dd=n-rl; 
wl=zeros (n) ; 

[Q  R]  =qr  (r)  ; 
e=eye (n) ; 
coun=0; 

for  i=n: -1 :n-rl+l; 

R11=R ( 1 : i , 1 : i) ;  %  FIND  THE  THE  LEADING  ixi  BLOCK  (STEP  1) 
[u, sigmin, v] =ssvd (Rll, 3 ) ;  %  FIND  THE  LEAST  DOMINANT 
% SINGULAR  VECTORS  AND  SINGULAR  VALUE  (STEP  2) 
pt=eye (i) ; 

vinf=norm(v, inf) ;  %  FIND  THE  MAXIMUM  ABSOL.  VALUE 

%  ELEMENT  OF  THE  LEAST  DOMINANT  SINGULAR  VECTOR 

vauxl=abs (v); 
ind=f ind (vauxl==vinf ) ; 

if  ind~=i,  %  FIND  THE  POSITION  OF  THE  MAX  ELEMENT 

vaux=pt  ( ind ,  :)  ;  %  MAKE  THE  PERMUTATION  TO  PUT  THE 

%  MAXIMUM  ELEMENT  AT  THE  i  th  POSITION  (STEP  3) 

pt (ind, : ) =pt (i, : ) ; 
pt (i, : ) =vaux; 

wl (:,i)=[v; zeros (n- i, 1) ] ;  %  ASSIGN  v  TO  THE  ith  COLUMN 

%0F  w  (STEP  4) 

p=pt ' ; 

ptil 12= zeros (i, (n-i) ) ; 
ptil21= zeros ( (n-i) , i) ; 
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%  STEP  5 
%  STEP  6 
%  STEP  7 


ptil22=eye( (n-i)  )  ; 
ptil=[p  ptill2  ;  ptil21  ptil22] ; 
w=ptil ' *wl; 

[Q1  Rlltilj  =qrgiv (Rll*p) ; 
e=e*ptil; 

R12=R ( 1 :  i ,  (i  +  1)  :n)  ; 
if  (n-i)==0, 

R2  2  =  [  ]  ; 
else 

R22=R ( (i+1) :n,  (i+1) :n) ; 

end; 

R= [Rlltil  Q1 ' *R12  ;  zeros ( (n-i) , i)  R22] ; 
Q12=zeros(i, (n-i) ) ; 

Q21=zeros ( (n-i) , i)  ; 

Q22=eye ( (n-i) ) ; 

Q=Q* [Q1  Q12  ;  Q21  Q22] ; 

end; 

end; 

nsbsp=Q( : ,n-rl+l:n) ; 
ssbsp=Q ( : , 1 :n-rl)  ; 


82 


Appendix  B 


function  jik  =  givensl (x,y,n, i,k) 

%GIVENS  Givens  rotation  matrix. 

%  G  =  GIVENS (x,y)  returns  the  complex  Givens  rotation 

matrix 
% 


%  j  c  s 

1  x 

r 

%  G  =  1 

such  that  G  * 

= 

%  J - conj ( s )  c 

0. 

I  y 

o 

where  c  is  real,  s  is  complex,  and  cA2  +  jsjA2  =  l. 


%  Copyright  (c)  1987-88  by  The  MathWorks,  Inc. 

%  Modified  by  Milton  P.  Ferreira  Sep/1992 

% 

%  Input  parameters : 

% 

%  x=  pivoting  element  [a(i,i)]. 

%  y=  element  we  want  to  zero  out  [a(k,i)] . 

%  n=  matrix  dimension. 

%  i=  column  of  the  element  we  want  to  zero  out. 

%  k=  row  of  the  element  we  want  to  zero  out. 

% 

%  Output  parameters: 

% 

%  j  ik=  matrix  "J".  When  multiplied  by  the  matrix  "a"  will  zero 
%  out  the  element  a(i,k) . 
absx  =  abs (x) ; 
if  absx  ==  0.0 

c  =  0.0;  s  =  1.0; 

else 

nrm  =  norm( [x  y] } ; 

c  =  absx/ nrm; 

s  =  x/absx*  (conj  (y) /nrm)  ; 

end 

j ik=eye (n)  ; 
jik(i, i) =c; 
jik(k, i) =-conj (s)  ; 
jik(i,k)  =s; 
j ik (k, k) =c; 
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Appendix  C 

%  Computes  Givens  Rotations  to  zero  out  the  elements  of  a 
%  matrix  below  its  main  diagonal .  It  is  used  at  step  6  of 
%  the  RRQR  algorithm  instead  of  a  new  QR  factorization  of 
%  Rll*P . 

%  Milton  P.  Ferreira  Sep/1992. 

% 

%  Input  parameters : 

% 

%  w=  matrix  whose  elements  below  the  main  diagonal  we  want  to 
%  zero  out. 

% 

%  Output  parameters: 

% 

%  ql=  new  matrix  Q,  after  applying  Givens  Rotations. 

%  rl=  new  matrix  R,  after  applying  Givens  Rotations. 

% 

function  [ql, rl] =qrgiv{w) 
i=size(w) ;i=i(lfl) ; 
qt=eye (i)  ; 
for  q=2:i; 

for  p=l:min( [q-1, i]  )  ; 
if  w (q, p)  ~=0, 

jpq=givensl (w(p,p) ,w{q,p) , i,p,q) ; 

qt=jpq*qt; 

w=jpq*w; 

end; 
end; 
end  ; 

ql=-qt ' ; 
rl=-w; 
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Appendix  D 


%  Algorithm  for  the  adaptive  tracking  of  a  moving  source 
%  Main  program. 

%  Milton  P.  Ferreira  Sep/1992. 

% 

%  Input  parameters: 

% 

%  db=  SNR. 

% 

%  Output  parameters: 

% 

%  R0=  autocorrelation  matrix. 

%  Rl=  noise- free  autocorrelation  matrix. 

%  Y=  each  column  of  Y  is  one  output  vector  "x"  from  the  "n" 

%  sensors. 

%  nest=  #  of  snapshots  used  to  compute  the  autocorrelation 
%  matrix. 

%  nupd=  #  of  updates  used  to  simulate  the  movement. 

%  ipp=  autocorrelation  matrix  dimension. 

%  nsin=  number  of  sources. 

% 

[RO , R1 , Y, nest , nupd, ipp,nsin] =cor4ml (db) ; 

[Q,  R, e, nsbsp, ssbspj =rrqr (Rl, ipp-nsin) ; 
for  i=nest+l :nest+nupd; 
old=Y ( : , i-nest)  ; 

[Q,R] =givqr (Q, R, e , old, -old) ; 
new=Y ( : , i )  ; 

[Q, R] =givqr (Q, R, e, new, new) ; 

[Q, R,  e, nsbsp,  ssbsp]  =rrqe  (Q,  R,  e,  nsin,  ipp)  ;  %  is  the  RRQR 
%  subroutine  without  the  initial  QR  decomposition 
[mags, angs] =minn (nsbsp, ssbsp, ipp) ; 

[ma, an] =ident (mags , angs) ; 
angupd ( i - nest ) =an (2,1) ; 

end; 

plot (i, angupd) ; 
grid; 

title ( ['SIGNAL  SUBSPACE  -  SNR= ' , int2str (db) , '  dB']); 

ylabel ( ' SOURCE  LOCATION' ) ; 

xlabel ( ' NB  OF  UPDATES ' ) ; 

text ( . 7 , . 5 , ' SOLID  =UPDATE ' , ' SC ' ) ; 
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Appendix  E 


%  Computes  the  rank-one  modification  of  a  square  matrix 
% 

%  Input  parameters: 

% 

%  Q=  orthonormal  matrix  resulting  from  the  RRQR 
%  factorization. 

%  R=  upper  triangular  matrix  resulting  from  the  RRQR 
%  factorization. 

%  e=  permutation  matrix  (PI) ,  resulting  fron  the  RRQR 
%  factorization. 

%  u  and  v=  u  and  v  vectors  that  modify  the  matrix  to  be 
%  updated. 

% 

%  Output  parameters: 

% 

%  Ql=  matrix  Q  after  modifying. 

%  Rl=  matrix  R  after  modifying. 

% 

%  Milton  P.  Ferreira  Sep.  1992. 

% 

function  [Q1,R1] =givqr (Q,R,e,u,v) ; 
w=Q ' *u ; 

Ql-Q; 

H=R; 

i=size(w) ; i=i (1,1) ; 
for  q-i-l: -1:1; 

jpq=givensl(w(i,l) ,w(q,l) ,i,i,q) ; 
w=jpq*w; 

H= jpq*H; 

Q1=Q1* jpq' ; 

end; 

Rl=H+w*v' *e; 
for  q«l:i-l; 

jpq=givensl(Rl(q,q) ,Rl(i,q) ,i,qfi) ; 

Rl=jpq*Rl; 

Ql=Ql*jpq'; 

end; 
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Appendix  F 


function  [mags,angs] =minn(nsbsp, ssbsp, ipp) 

%  compute  the  noise  and  signal  zero  locations  using  the  TK 
%  min -norm 

%  alg.  (version  1.0  10/14/91),  Monique  P.  Fargues. 

% 

%  Input  parameters: 

% 

%  ipp=  correlation  matrix  dimension 
%  nsbsp=  noise  subspace. 

%  ssbsp=  signal  subspace. 

% 

%  Output  parameters: 

% 

%  mags=  magnitude  of  the  polynomial  roots . 

%  angs=  phase  of  the  polynomial  roots. 

% 

% 

dn=zeros (l:ipp) ;ds=zeros (l:ipp) ; 
clear  g 

g=nsbsp ( 1 , : ) ;  %noise  zeros 

En=nsbsp (2 : ipp, : ) ; 

dn(2:ipp) =En*g' / (g*g' ) ; 

dn (1) =1; 

clear  g 

g=ssbsp (1, : ) ;  %signal  zeros 

Es=ssbsp (2 : ipp, : ) ; 

Kg=- 1/ (l-g*g' ) ; 

ds (2 : ipp) =Kg* (Es*g' ) ; 

ds (1) =1; 

flops (0) ; 

dnr= roots (dn) ; 

dsr= roots (ds) ; 

mags=abs (dsr) ; 

angs=angle (dsr) *180/pi ; 
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Appendix  G 


%  Computes  the  "refined"  signal  subspace 
%  Version  1.0  08/21/92,  Monique  P.  Fargues. 

% 

%  Input  parameters 
% 

%  R=  upper  triangular  matrix  resulting  from  the  RRQR 
%  factorization. 

%  e=  permutation  matrix  (PI)  resulting  from  the  RRQR 
%  factorization. 

%  nsin=  #  of  sources. 

%  ipp=  correlation  matrix  dimension. 

% 

%  Output  parameter: 

% 

%  ref=  "refined"  and  reorthonormalized  signal  subspace. 
% 

function  ref =ref ine (R, e, nsin, ipp) 
pil=e ( : , lrnsin) ; 
pi2=e ( : , nsin+1 : ipp) ; 

Rll=R(l:nsin, l:nsin) ; 

R12=R(l:nsin, nsin+1: ipp) ; 
ref=orth (pil*Rll ' +pi2*R12 ' )  ; 


88 


Appendix  H 


%  identify  the  signal  magnitude  and  location  of  the  roots 
%  version  1.0  10/13/91,  Monique  P.  Fargues. 

%  Input  parameters : 

% 

%  mag=  magnitude  of  the  polynomial  roots 
%  ang=  phase  of  the  polynomial  roots 
% 

%  Output  parameters: 

% 

%  mag_r=  vector  containing  the  magnitude  of  the  roots 
%  corresponding  to  the  sources . 

%  ang_r=  idem  to  the  phases. 

% 

function  [mag_r, ang_r] =ident (mag,ang) 

[m,n] =size (ang) ; 
k=l; 

ang_r0  (1:2 , 1)  *  [0;  0]  ;mag_r0  (1 : 2 , 1)  =  [0 ;  0]  ; 
for  i=l:m 

if  mag(i)<1.3  &  ang(i)<50  &  ang (i) >10  %look  at  both 
ang_r0 (k, 1) =ang(i)  ; 
mag_r0 (k, 1) =mag  ( i ) ; 
k=k+l ; 
end 

if  k==3 ,  break,  end 
end 

%  sort  to  insure  proper  separation  of  sines 
[ang_r ( : , 1) , 1] =sort (ang_r0 ( : , 1) ) ; 
mag_r (:,!)= (mag_r0 (1,1)); 
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function  [sigmin,x] =icest (R) 

%  1/7/92  ******-icest .m-****** 

%  version  1.0,  Monique  P.  Fargues. 

%  incremental  condition  estimator  (modified  from  Bischoff 
%paper) 

%  computes  an  estimate  for  minimum  singular  value  and 
%vector 
%  associated 

%  with  an  upper  triangular  matrix 
%  initial  matrix  gamma  v' 

%  OR 

%  Input  parameters : 

% 

%  R=  matrix  we  want  to  estimate  the  smallest  S.  value  and 
%  its  S.  vectors. 

% 

%  Output  parameters: 

% 

%  sigmin:  minimum  singular  value 

%  x:  .  singular  vector 

%  . - . - . 

clear  xx  x  v 

[m,  n]  =size  (R)  ; 

%if (m~=n) ,  error  CR  is  not  square'},  end 

%if (any (diag (R) ==0) )  %  matrix  is  singular 

%  smin=0;  vmin=zeros (n, 1) ;vmin (min (find (diag (R) ==0) )) =1; 

%  return 

%end 

% 

x(n, 1) =1/R(n,n) ; 

for  i=n-l:-l:l 

xx(l:n-i,l)=x(i+l:n,l) ; 
v=R(i,i+l:n) ' ;  gamma=R ( i , i ) ; 
alpha=v' *xx; 
if  alpha~=0 

beta= (abs (gamma) *norm(xx, 2) ) a2  +  abs (alpha) A2  -1; 
eta=beta/ (2*abs (alpha) ) ; 
sqr=sqrt (etaA2+l)  ; 

nu=eta  +  sqr  ;  %  sqrt  (eta>s2+l)  ; 

temp=abs (alpha) *nu; 
root=temp+l;  sqr=sqrt (nuA2+l) ; 

s=temp/ (alpha*sqr)  ;  %  sqrt (nuA2+l) ) ; 

c=-l/sqr  ;  %  sqrt  (nuA‘2+l)  ; 


else 


root= (abs (gamma) *norm(xx,2) ) *2; 
if(root>l),  c=l ,  s=0 ; 
else,  c=0,  s=l ;  end 
end 

x(i:n, 1) = [ (c-s*alpha) /gamma; s*x (i+1 :n, 1) ] ; 

end  %  of  initial  loop 
xnorm=norm(x(i :n, 1) ,2) ; 
x(l:n, 1) =x(l:n, 1) /xnorm; 
s igmin= 1 / xnorm ; 
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Appendix  .1 


function  [vsv, sigmin, usv] =ssvd(a,k) 

%  THIS  FUNCTION  IS  THE  IMPLEMENTATION  OF  THE  INVERSE 
% ITERATION 

%  TECHNIQUE  TO  FIND  THE  SINGULAR  VECTORS  AS  SHOWN  IN  THE 
% PAPER 

%  "DEFLATED  DECOMPOSITION  OF  SOLUTIONS  OF  NEARLY  SINGULAR 
%  SYSTEMS"  -  TONY  F.  CHAN  -  PG  746  -  SIAM  J.  NUMER.  ANAL. 
%  VOL  21  #4  AUG.  1984 
%  [vsv, sigmin, usv] =ssvd (a, k) . 

% 

%  Input  parameters : 

% 

%  a=  matrix  we  are  looking  for  the  smallest  singular 
%  vectors. 

%  k=  number  of  iterations  desired. 

% 

%  Output  parameters: 

% 

%  usv,  vsv=  right  and  left  smallest  S.  vectors. 

%  sigmin=  minimum  singular  value. 

% 

%  Milton  P.  Ferreira 
% 

n=size (a) ; 
n=n(l, 1) ; 
v=ones (n, l) ; 
for  i=l:k; 

util=a\v; 

u=util/norm(util) ; 
vtil=a' \u; 
v=vtil/norm(vtil) ; 

end; 

vsv=v; 

util=a\vsv; 

usv=util/norm(util) ; 

sigmin=l /norm (util) ; 
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Appendix  K 


function  [R0,R1, Y,nest,nupd, ipp,nsin] =cor4ml (db) ; 

% 

%  Monique  P.  Fargues. 

% 

%  COMPUTE  THE  CORRELATION  FUNCTION  ONLY 
%  RO:  Toeplitz  correlation  function 
% 

%  Input  parameters : 

% 

%  db=  SNR. 

% 

%  Output  parameters: 

% 

%  R0=  autocorrelation  matrix. 

%  Rl=  noise- free  autocorrelation  matrix. 

%  Y=  each  column  of  Y  is  one  output  vector  "x"  from  the  "n" 
%  sensors. 

%  nest=  #  of  snapshots  used  to  compute  the  autocorrelation 
%  matrix. 

%  nupd=  #  of  updates  used  to  simulate  the  movement. 

%  ipp=  autocorrelation  matrix  dimension. 

%  nsin=  number  of  sources. 


clg; format  compact 
seedl=1042 ; rand ( ' seed' , seedl) 

icor=0;  %input ( ' true/est  correl  1/0:  '); 
itop=0;  %input ( ' toeplitz/non  toeplitz  1/0: 

nupd=200;  %input ( 'update  nb  nupd:  '); 
%nupdO=input ( 'update  nb  nupdO  (drop  in  angle):  '); 
del_freq=5;  %  input (' del_freq  (in  percen*f req (1) : 
%if  icor==l 
%  fprintf ( ' 

%else 

%  fprinrf ( ' 

%end 

nest=100 ; 

%gdb= input ( ' 
gdb= [db, db] ; 
ang=[30  40]; 
freqt=[7.5  9] 
ipp=10; 
nsin=2 ; 
ssig=l . ; 
jc=sqrt ( 
sigma=l ; 


true  cor.  seq\n' 
cor.  seq\n' 


est 


input  gdb:  [x  x]  ') 


% 

% 

% 

% 

% 


-) 


source  angles 
temporal  frequencies 
number  of  sensors 
number  of  sources 
ref.  noise  variance 
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%  linear  amplitude 


for  i=l:nsin 

A(i) =sqrt (2 . ) * (10A (gdb (i) /20 . ) ) ; 
end 

if  icor==l 

%  compute  the  true  correlation  sequence 
for  k=l:ipp 
res=0; 
if  k— l 
res=ssigA2 ; 
end 

for  i=l:nsin 

res=res+ ( A { i ) A2) *exp ( jc* (k-1) *ang (i) *pi/180 . ) /2 . ; 
end 

R (k) =res; 
end 

RO=toeplitz (R) ; 
else 

%  compute  the  estimated  correlation  seq 
%  based  on  nest  data  points 
rand ( ' normal ' ) 

%sigm=sqrt (10A ( (sigma) /10) ) ; 
for  i=l:nsin 

A(i) =sqrt (2) *sigma* (10A (gdb (i) /20) ) ; 
end 

for  i=l:ipp 

Y (i, 1 :nest+nupd) =sigma* (rand (1, nest +nupd) +jc*rand (l.nest+nup 

d)  )  ; 

end 

f req=ang*pi/180 ; 

rand (' uniform' )  %  create  uniform  variable  dist. 

%  in  (-pi, pi) 

Xl=pi* (rand(l,nsin*ipp* (nest+nupd) )-.5); 
freqO=freq (1) ; 

for  j 2=1 : nest+nupd  %  j 2 :  time 

snapshot 

%f req (1) =f reqO  -  (pi/180) * (j2-l) *del_f req/nupd;  %del_freq 

%  change  in  nupd  samples 

%if  j2>nest+nupd0 ,  f req (1) =f reqO*del_f req;  end  %step  freq. 

%  change 

for  i=l:nsin  %  i:  number  of  sines 

for  jl=l:ipp  %  jl:  sensor  position 

temp= ( jl*freq(i) +j2*freqt (i) +X1 (1, j2+ (i-1) *ipp) ) ; 

Y ( j 1 , j  2 ) = Y ( j 1 , j  2 ) +A ( i ) *exp ( jc*temp) ; 
end 
end 
end 

R0=Y ( : , l:nest) *Y ( : , l:nest) ' ; 
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end 


%get  the  noise  free  version  of  RT  &  normalize  the  matrix 
%  compute  the  EVD  of  the  noisy  matrix  for  later  computations 


%  for  original  correlation 
% [Ve,De] =eig (RO)  ; 

%  sort  the  eigenvalues 
% [Des , 1] =sort (diag (De) ) ;12 
%Ves=Ve ( : , 12 ( : ) ) ; 

%Vnoi=Ves { : , nsin+1 : ipp) ; 
%Vsig=Ves ( : , 1 :nsin)  ; 
Rl=R0-nest*ssigA2*eye (RO) ; 
matrix 

%R1=R1 . /R1 (1,1)  ; 


function 


flipud (1)  ; 

%  sorting  in  descending  order 
%  isolate  the  noise  vectors 

%  . -  signal  . . 

%  noise  free  correlation 

%  potential  matrix 
%  normalization 


if  icor==l, 

Rl=R0-ssigA2*eye (RO) ;  %  noise  free  correlation 

matrix 

end; 

if  icor==0 , 

Rl=R0-nest*ssig>k2*eye  (RO)  ; 

end; 

if  itop==l , 

Rl=toep (Rl) ; 

end 
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Jlo  *XO  -V?  Jlo  nV?  «iO  Jka  An  An  -Wi  jm  An  An 
©V*  Or  OT"  Or  ©v'  0V'  Ov'  ©V'  Or  ©’r'  ©T'  «r  Crr' 


Transforms  a  matrix  to  a  Toeplitz  matrix.  Used  with 
"cor4ml .m" . 

Input  parameter: 

Matrix  to  be  transformed. 

Output  parameter: 

sum=  Toeplitz  matrix  after  transformation. 

Milton  P.  Ferreira  Sep/1992. 

function  sum-toep (Rl) ; 
n=size(Rl) ;n=n(l,l) ; 
sum= zeros  (Rl)  ; 
for  i=-n+l:n-l; 

a=diag (Rl, i) ; 
s=mean(a) ; 
sl=size (a) ; 
sl-sl(l,l) ;  . 
s2=ones (si, 1) *s; 
sum=sum+diag (s2 ,  i)  ; 


end; 


<#><*>(#><#><#><#><#>!*><#><#><#>«*><#>'#> 


Appendix  M 


Computes  the  smallest  eigenvector  by  inverse  iteration. 
Input  parameters: 

a=  matrix  from  which  we  want  to  compute  an  approximation 
of  the  smallest  eigenvector. 
u=  #  of  iterations. 

Output  parameter: 

x=  approximation  for  the  smallest  eigenvector. 

Milton  P.  Ferreira  Sep/1992. 

function  x=eeig(a,u) 
n=size (a) ; 
n=n (1, 1)  ; 
x=ones (n, 1)  ; 
for  i=l:u; 
t=a\x; 

x=t/norm(t , inf ) ; 

end; 

x=x/norm(x)  ; 
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